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1.0  OBJECTIVE 


The  research  performed  under  this  contract  is  directed  toward  an  under- 
standing of  plasma  processes  in  boundary  and  electrode  regions.  The 
aim  is  to  develop  a self-consistent  theory  combining  electron  and  ion  kinet- 
ics in  the  presence  of  strong  field  and  density  gradients  with  the  boundary 
layer  equations  of  gasdynamics.  The  present  study  has  direct  application 
to  problems  in  high-pressure,  flowing  gas  discharges  and  in  MHD  power 
generators . 

In  electric  discharges,  the  development  of  arcs  and  nonuniform  current 
distributions  begins  at  the  electrodes.  The  interaction  between  a thermal 
boundary  layer  and  the  space-charge  sheath  adjacent  to  electrodes  and  other 
solid  boundaries  is  therefore  critical  to  the  determination  of  electric  discharge 
stability  and  uniformity.  If  an  understanding  of  the  basic  processes  in  these 
regions  can  be  achieved,  it  may  be  possible  to  extend  the  useful  operating 
range  of  discharge-pumped  devices  to  higher  power  densities  and  longer 
pulses . 

In  an  MHD  generator,  the  internal  voltage  losses  across  the  boundary 
layer  and  sheath  limit  the  energy  conversion  efficiency.  Current  con- 
strictions on  the  electrodes  can  cause  surface  deterioration.  An  under- 
standing of  the  collisional  sheath  and  boundary  layer  problems  will  lead 
to  improved  system  designs. 
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2. 0 SUMMARY 

In  the  first  quarter,  a literature  search  was  performed  to  identify  relevant 
theories  and  to  evaluate  their  applicability  to  the  present  study.  The  important 
processes  which  must  be  included  in  a self-consistent  model  of  the  plasma 
sheath  were  identified.  The  Boltzmann  equation  for  electrons  in  a nonuniform 
field  was  derived  in  a form  amenable  to  numerical  solution. 

In  the  second  quarter,  a computer  program  was  developed  to  model  the  electron 
kinetics  in  a highly  nonuniform  electric  field.  Cases  were  run  in  several 
| gases  using  a constant  field  for  comparison  with  previous  Monte  Carlo  solu- 

tions. The  nonequilibrium  character  of  the  electron  energy  distribution  was 
clearly  demonstrated,  thus  pointing  out  the  limitations  of  earlier  sheath 
models.  The  equations  for  electric  field  and  ion  kinetics  were  formulated  for 
use  in  a fully  self-consistent  solution  of  the  cathode  sheath. 

In  the  third  quarter,  the  evolution  of  electron  energy  distributions  in  the 
cathode  fall  and  negative  glow  was  examined  by  assuming  a linearly  falling 
electric  field  followed  by  a zero-field  drift  region.  The  results  agree  well 
with  experimental  distributions  measured  at  the  edge  of  the  cathode  fall. 
Recombination  and  attachment  were  included  in  the  model  but  only  had  a sig- 
nificant effect  in  the  negative  glow.  Analytic  solutions  to  the  Boltzmann 
equation  were  derived  for  a specific  choice  of  cross  sections  in  order  to 
test  some  of  the  angular  approximations  used  in  the  code. 

In  the  fourth  quarter,  cathode  sheath  characteristics  were  derived  using  the 
linear  field  approximation.  It  was  found  that  the  usual  scaling  laws  apply 


in  terms  of  the  reduced  pressure  at  the  cathode  surface.  Gas  heating  in  the 
sheath  was  examined  and  the  temperature  distribution  was  calculated.  For 
subsonic  flows,  the  principal  cooling  mechanism  is  conduction  to  the  cathode 


surface.  In  this  regime,  the  boundary  layer  problem  can  be  separated  from 
the  sheath  problem.  A computer  program  for  the  solution  of  the  boundary  layer 
equations  including  gas  heating  was  described. 

In  the  last  quarter,  numerical  calculations  of  the  temperature  distribution 
in  the  cathode  sheath  and  boundary  layer  were  presented.  The  temperature 
rise  across  the  sheath  was  shown  to  increase  linearly  with  pressure.  An 
analysis  was  made  of  the  thermal  stability  at  the  cathode  surface.  A critical 
current  density  was  found  at  which  the  surface  temperature  and  electric 
field  increase  until  arcing  occurs.  This  event  can  be  forstalled  by  cooling 
the  electrode  and  by  the  use  of  a buffer  gas  with  high  thermal  conductivity. 
Auxiliary  ionization  and  field  tailoring,  which  are  designed  to  eliminate  ion- 
ization instabilities  and  maintain  constant  E/N  in  the  positive  column,  affect 
the  sheath  only  through  the  imposed  local  current  density  and  do  not  sig- 
nificantly alter  sheath  characteristics. 


3.  0 


INTRODUCTION 


The  problem  as  stated  in  Section  1.0  is  to  develop  a self-consistent  theory 
of  plasma  processes  in  boundary  and  electrode  regions.  The  approach  taken 
under  this  contract  consists  of  two  basic  parts:  First,  a refinement  of  the 
existing  theories  of  the  sheath  and  cathode  fall  regions  to  include  the  effects 
of  field  and  density  gradients  and  the  anisotropy  of  the  electron  velocity  dis- 
tribution; then,  a coupling  to  the  gas  dynamic  equations  with  transverse  flow 
to  study  the  growth  of  boundary  layer  thickness  and  the  effect  of  convective 
cooling  on  sheath  structure.  The  basic  processes  which  must  be  included  in 
this  model  are  described  below. 

A plasma  consists  of  electrons,  ions,  and  neutral  particles  in  the  presence 
of  electric  and/or  magnetic  fields.  If  the  gas  is  weakly  ionized,  and  all 
collisions  are  binary  encounters,  the  particles  are  completely  described  by 
distribution  functions  which  satisfy  a Boltzmann  equation  with  the  appropriate 
Lorentz  force  and  collision  term.  In  most  practical  cases,  the  distribution 
function  for  neutral  particles  is  close  to  the  equilibrium-Maxwellian  form, 
so  that  moment  equations  can  be  used  to  describe  the  local  gas  properties 
in  terms  of  an  average  density,  flow  velocity  and  temperature.  For  a sheath 
imbedded  in  a turbulent  boundary  layer,  the  effects  of  viscosity  and  thermal 
conductivity  must  be  included  along  with  a source  term  for  gas  heating. 

Because  of  their  mass,  the  ions  can  also  be  described  by  moment  equations, 
with  the  inclusion  of  relevant  production  and  loss  processes,  momentum  and 
energy  transfer  to  neutrals,  and  Joule  heating.  However,  due  to  their  greater 
mobility  and  resulting  nonequilibrium  behavior,  the  electrons  must  be  dealt 
with  in  terms  of  the  full  space  and  velocity  distribution  function.  Electrons 
play  a dominant  role  in  determining  sheath  structure.  Consequently,  their 
accurate  description  is  considered  an  essential  part  of  the  present  study. 

The  electron  velocity  distribution  must  be  obtained  from  a numerical  solution 
of  the  colli sional  Boltzmann  equation  without  the  usual  angular  approximation. 
The  treatment  of  electron  and  ion  kinetics  must  be  self-consistent  with  the 
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solution  of  Poisson's  equation  for  the  electric  field.  In  the  presence  of  an 
applied  magnetic  field,  additional  field  equations  are  required  along  with  a 
modified  Lorentz  force  term. 


3.  1 Literature  Search 

A literature  search  was  conducted  covering  the  areas  of  nonequilibrium 
electron  kinetics  and  collision-dominated  sheaths.  The  list  of  references 
presented  is  meant  to  be  more  representative  than  exhaustive.  Preference 
was  given  to  the  most  recent  papers  and  the  latest  work  of  a given  author. 
Studies  of  the  turbulent  boundary  layer  which  did  not  treat  electric  sheaths 
were  not  included.  The  considerable  research  on  electric  probe  theory 
was  also  not  considered  immediately  applicable.  Only  one  paper  was 
found ^ which  analyzed  the  interaction  between  sheath  and  boundary  layer. 
The  results  of  that  work  will  be  incorporated  into  the  present  study. 


The  term  "nonequilibrium  electron  kinetics"  refers  to  either  spatial  or 
temporal  nonequilibrium  as  occurs,  for  example,  in  the  Townsend  dis- 
charge and  in  streamer  development,  respectively.  In  either  case,  the 
electron  energy  distribution  is  not  completely  specified  by  the  local  elec- 
tric field  and  gas  density  but  must  be  considered  in  the  full  context  of 

spatial  or  temporal  growth.  Studies  of  electron  kinetics  can  be  broken 

, ...  5,  10,  12,  13,  15 

into  two  major  groups,  1)  those  based  on  Monte  Carlo  simulation  . 

. 4,  6, 7,  11,  14,  16 

and  2)  those  based  on  analysis  of  the  Boltzmann  equation. 

4,  5,  8 

They  include  several  detailed  comparisons  of  the  two  methods.  ’ ’ Most  of 

the  calculations  were  based  on  a uniform  electric  field,  except  for  one  assuming 

9 10 

a logarithmic  and  one  a linear  dependence.  There  were  only  two  papers 

involving  the  self-consistent  treatment  of  nonequilibrium  electron  kinetics 

with  the  field  equations.  One  was  a Monte  Carlo  simulation  of  the  temporal 

13 

development  of  electron  avalanches  or  streamers.  The  other  was  a study 

of  a low-pressure  plasma  diode  using  the  Krook  model  to  describe  the  col- 

14 

lision  term  in  the  Boltzmann  equation.  These  theories  will  be  discussed  in 
the  following  section. 


, 17 


The  subject  of  collision-dominated  sheaths  was  broadened  to  include  other 
situations  involving  the  interaction  between  a plasma  and  space-charge- 

induced  fields.  The  papers  in  this  category  can  be  divided  roughly  into 

, , u lo.  20, 24, 26-30 

two  groups,  1)  analyses  of  steady-state  sheaths 


and  2)  simula- 


tions of  discharge  initiation  ’ and  glow-to-arc  transitions.  Both  prob- 

lems involve  the  same  system  of  equations  and  there  is  little  difference  in  the 
general  method  of  solution.  Only  the  assumptions  made  and  the  proces  ses  included 
differentiate  these  studies.  All  treatments  are  based  on  the  use  of  moment 

equations  to  describe  the  electrons  and  ions.  Some  consider  gas  heating 

21  , 20,  24, 27  , . 

effects  . Others  include  an  external  ionization  source.  Most  models 

deal  with  only  one  spatial  dimension  (except  References  18,  21,  and  23, which 

consider  two).  The  applicability  of  these  approaches  to  the  present  study  will 

be  addressed  in  the  following  section. 


3 . 2 Evaluation  of  Appropriate  Theories 

The  theories  identified  in  the  literature  search  and  listed  in  Section  3.4 
will  now  be  evaluated  on  the  basis  of  their  potential  in  meeting  the  require- 
ments of  the  present  study.  The  requirements  are  as  follows: 

a)  The  system  of  equations  must  include  as  a minimum,  the 
Poisson  equation  for  the  electric  field,  the  collisions!  Boltz- 
mann equation  for  electrons,  the  moment  equations  for  ions 
and  the  boundary  layer  equations  for  the  neutral  gas.  The 
kinetic  equations  for  electrons  and  ions  may  be  replaced  by 
Monte  Carlo  simulation  with  an  appropriate  collision  model. 

b)  The  theory  must  consider  the  proper  boundary  conditions  at 
the  wall-sheath  interface  and  at  the  sheath-plasma  interface. 

The  sheath-plasma  boundary  conditions  may  be  omitted 
if  the  calculations  are  continued  to  the  opposite  wall  or  elec- 
trode. 
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c)  The  solution  must  be  obtainable  in  two  spatial  dimensions 
with  full,  self-consistent  coupling  between  all  equations. 

d)  The  treatment  of  electron  kinetics  must  utilize  realistic 
cross  sections,  including  angular  dependence,  a source 
term  for  external  ionization  and  loss  terms  for  attach- 
ment and  recombination. 

e)  The  equations  for  ion  kinetics  must  include  the  inertia 
term,  charge  exchange,  and  the  same  production-loss 
terms  as  for  electrons. 

f)  The  boundary  layer  equations  must  consider  gas  heating 
in  the  sheath,  as  well  as  heat  conduction  to  the  wall. 

g)  Since  the  computer  code  is  to  be  used  in  a parametric 
mapping  of  sheath  characteristics  over  quite  a broad 
parameter  space,  it  must  be  reasonably  economical  in 
the  use  of  central  processor  memory  and  time. 

First,  a comparison  will  be  made  between  Boltzmann  and  Monte  Carlo 
techniques.  Although  the  two  approaches  are  mathematically  equivalent, 
there  are  important  economic  and  scientific  differences.  Economically, 
the  Monte  Carlo  methods  are  very  costly  as  far  as  computer  time  is  con- 
cerned. No  quantitative  comparisons  were  found,  but  all  authors  who 
tried  both,  found  the  Monte  Carlo  simulation  to  be  much  more  time  con- 
suming than  a Boltzmann  analysis.  Scientifically,  the  Boltzmann  approach 
gives  a greater  insight  into  the  physics  of  a particlar  situation.  Approx- 
imate solutions  can  be  found  analytically  which  illustrate  the  dependence 
on  relevant  parameters,  and  the  sensitivity  to  changes  in  input  data  can 
be  more  easily  determined.  Also,  exact  analytic  solutions  in  certain 
situations  provide  a benchmark  for  numerical  results.  The  principal 
advantage  of  the  Monte  Carlo  method  is  the  ease  with  which  it  can  be  form- 
ulated into  a computer  code. 
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Recent  solutions  of  the  Boltzmann  equation  have  dealt  with  the  problem  of 

electron  swarm  motion  in  a uniform  electric  field,  when  ionization  and/or 

attachment  are  appreciable.4’  ^ These  studies  treat  the  effect  of  electron 

density  gradients  by  an  expansion  in  higher-order  spatial  derivatives. 

(The  present  author  has  calculated  the  diffusion  coefficients  for  these 

higher-order  terms  in  argon  at  low  values  of  E/N,  where  electron  impact 

3 6 

ionization  is  negligible.  ) Using  this  approach,  the  electron  density  in 
an  avalanche  or  swarm  can  be  described  by  an  expanded  continuity  equa- 
tion with  position-  and  time-  independent  coefficients.  A comparison  with 

5,  8 

Monte  Carlo  results  indicates  excellent  agreement. 


The  technique  described  above  is  not  immediately  applicable  to  the  sheath 
problem,  because  of  the  assumption  of  a uniform  electric  field.  An  exten- 
sion of  the  method  to  include  spatial  gradients  of  the  field  does  not  appear 

, 7 

practical.  A more  feasible  approach,  outlined  by  Allis,  involves  a 
direct  solution  of  the  Boltzmann  equation  without  expansion.  He  derives 
an  equation  for  the  flux  of  electrons  in  the  two-dimensional  space  (W,  V), 
where  W is  the  total  electron  energy  and  V is  the  electric  potential.  This 
approach  neglects  the  sjecific  angular  dependence  of  the  distribution  func- 
tion and  hence  prevents  the  use  of  differential  cross  sections.  We 
propose  to  include  a third  independent  variable  defining  the  component  of 
electron  motion  in  the  field  direction  and  to  solve  the  Boltzmann  equation 
numerically.  This  method  will  be  discussed  in  the  following  section. 

A technique  similar  to  the  one  proposed  here  has  been  used  by  Fournier, 
et  al^  to  analyze  the  cathode  region  of  a glow  discharge.  His  work  has 
not  been  published  in  detail,  so  it  is  not  clear  what  calculational  methods 
were  used  or  what  assumptions  were  made.  Preliminary  results 
indicate  that  this  method  provides  an  accurate  description  of  the  electron 

3 9 

kinetics  in  a sheath.  Other  approaches  ’ which  make  simplifying  assump- 
tions in  order  to  ease  the  numerical  calculation  are  useful  in  understanding 
the  physics  of  a situation  and  in  arriving  at  analytic  expressions  for  trans- 

8 


4 

A 


port  and  ionization  coefficients  but  are  not  considered  accurate  enough  or 
general  enough  for  the  present  study. 

An  analysis  which  comes  very  close  to  meeting  the  requirements  of  our 

14 

model  for  electron  and  ion  kinetics  was  presented  by  Whipple,  et  al. 

They  obtained  a self-cons' stent  solutionof  the  Boltzmann  equation  for 
electrons  and  ions  along  with  Poisson's  equation,  in  a low-pressure  plasma 
diode.  The  main  shortcoming  of  this  work  was  the  use  of  an  inverse  fifth 
power  law  potential  to  describe  collisions  between  charged  particles  and 
neutrals.  The  method  of  solution,  however,  involves  a nested  iteration 
scheme  which  should  work  equally  well  with  the  more  accurate  collision 
term  required  here. 

The  analyses  of  collision-dominated  sheaths  cover  a variety  of  special 
situations,  from  MHD  generators  to  gas  discharges,  from  the  cathode 
fall  region  to  the  positive  column,  both  transient  and  steady  state.  The 
equations  used  to  describe  the  plasma,  however,  differ  little  from  one 
treatment  to  the  next.  The  methods  of  solution  depend  on  geometry  and 
time  scale,  but  qualitatively  are  very  similar.  The  principal  drawback 
in  all  cases  is  in  the  treatment  of  electron  kinetics.  The  use  of  moment 
equations  with  equilibrium  coefficients  is  more  accurate  in  some  situations 
than  in  others  but  always  falls  short  of  providing  a complete  model. 


What  all  the  studies  have  shown  is  that  in  treating  sheaths  and  streamers, 
it  is  important  to  consider  neutral  density  variations,  space  charge  induced 
fields  and  nonequilibrium  electron  kinetics.  However,  these  effects  have 
never  been  brought  together  into  a single,  unified  theory  of  the  plasma  sheath. 

The  only  significant  addition  to  the  model  of  the  cathode  fall  since  Ward's 

, , 29,  30  . . . * , . , , 20,  24, 27 

early  work  has  been  the  inclusion  of  an  external  source  term. 

2 3 6 

The  recent  renewal  of  interest  in  the  problem  ’ ’ should  lead  to  signifi- 

cant contributions  in  the  near  future. 


MATHEMATICAL  MODEL 


4.  0 

The  equations  describing  electron  and  ion  kinetics  in  an  electric  field  are  pre- 
sented below.  The  Boltzmann  collision  term  for  electrons  is  derived  in  its 
full  angular  dependent  form  so  as  to  be  applicable  in  the  case  of  strongly 
anisotropic  distributions.  Elastic  and  inelastic  collisions  are  treated  along 
with  various  electron  production  and  loss  processes.  The  continuity  and 
momentum  equations  for  positive  and  negative  ions  are  modified  to  include  a 
source  term.  The  gas  dynamic  equations  will  be  presented  in  Section  5.  2. 

4.  1 The  Boltzmann  Equation 

The  distribution  function,  f(c7  "?T  t ).  provides  a complete  statistical  descrip- 
tion of  the  electrons  in  a weakly-ionized  gas.  The  number  of  electrons  in 
the  volume  element  (r,  r+dr)  with  velocities  in  the  range  (?T  c+dcT)  at  time 
t is  given  by  f(c,  r,  t)  dc  dr  dt.  The  equation  which  describes  the  evolution 
of  this  function  in  space  and  time  is  the  collisional  Boltzmann  equation, 

TT  + ~ V + r ' V ’ (-If)  <» 

c 

where  a is  the  acceleration  due  to  electric  and  magnetic  fields.  In  general, 

a”  = - — ( E + c x B)  (2) 

m 

where  e is  the  absolute  value  of  the  electronic  charge.  The  term  on  the 
right  in  Equation  (1)  represents  the  change  in  f due  to  collisions  and  will  be 
described  in  detail  below. 

If  we  consider  for  the  moment,  variations  in  one  spatial  dimension,  with 
axis  antiparallel  to  the  applied  electric  field,  then  Equation  (1)  can  be 
rewritten  for  f(c,  c^,  y)  as 
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where  E(y)  is  the  magnitude  of  the  local  electric  field.  By  dividing  both  sides 


of  Equation  (3)  by  c^  and  eE(y)  we  have, 


di  df 

<30  + d£ 


eCyE(y) 


(SB 


(4) 


where 


r y 1 2 

0 = e/  E(y)dy  is  the  electron  potential  energy  and  £ = y mc^ 
is  the  electron0  kinetic  energy  in  the  field  direction.  If  we  define  the 
total  kinetic  energy  as  ( = y mc^,  then  f can  be  represented  as  the  sum 
of  two  functions,  f (€,£,0)  and  f (6,4,0),  for  the  distribution  of  electrons 
moving  with  ( cy  < 0)  and  against  ( Cy  > 0)  the  electric  field,  respectively. 

In  the  absence  of  collisions,  f is  constant  along  the  lines  4 ■ 0 = const. 

These  are  known  as  the  characteristics  and  also  correspond  to  electron 
trajectories  in  the  two-dimensional  space  (4,0).  When  collisions  are 
included,  the  term  on  the  right  in  Equation  (4)  represents  the  change  in 
f along  the  characteristics.  As  we  shall  see  below,  the  collision  term  is 
a complicated  function  of  f and  the  various  cross  sections.  The  present 
method  of  solution  will  be  to  estimate  f,  use  this  in  evaluating  the  colli- 
sion term,  and  then  integrate  Equation  (4)  along  the  characteristics.  The 
distribution  function  thus  obtained  becomes  the  next  estimate  in  an  iterative 

procedure  which  continues  until  convergence  is  reached.  This  technique 

31  3334 

has  been  used  in  similar  situations  by  the  present  author  and  others 
with  considerable  success. 


The  boundary  conditions  for  f are  given  at  0 = 0 and  those  for  f at 
0=0  . The  function  f,(6,  4,0)  may  represent  the  distribution  of 

electrons  ejected  from  a cathode  surface  by  ion  impact,  or  it  may  equal 
zero  for  a nonemitting  surface.  The  function  f (€,  4,  0 ) is  either 


max 


the  distribution  of  electrons  at  the  head  of  the  positive  column  (0 


max 


or  at  the  anode  ( f_  = 0),  or  along  an  axis  of  symmetry  (f  = f+).  Along 
the  locus  of  turning  points,  4 = 0,  the  condition  f ( €,  0, 0 ) =f  (6,0,0) 
applies . 


II 


Boltzmann  Collision  Term 


! 


/ 


I 


I 


The  collision  term  in  the  Boltzmann  equation  represents  the  rate  at  which 
electrons  are  scattered  into  or  out  of  a given  region  of  phase  space  as  a 
result  of  collisions  with  ions,  neutrals,  or  other  electrons.  For  the  pres- 
ent, we  shall  consider  only  collisions  with  neutral  p4- tides  and  recom- 
bination with  positive  ions.  The  neutral  particle  reactions  include  elastic, 
inelastic  and  superelastic  scattering,  ionization  and  attachment.  The 
first  three  processes  result  in  a change  of  phase  space  coordinates  for  a 
single  electron,  whereas  the  second  two  processes  involve  the  creation  or 
removal  of  an  electron.  Positive  ion  recombination  also  involves  the  loss 
of  an  electron.  The  collision  term  is  then  the  sum  of  contributions  from 
the  various  processes,  i.  e. 


Ha  ■ m • (*)  * m 


(jl\ 

\at/ 


(-if) 

a,  r 


(5) 


The  individual  terms  will  now  be  derived, 
a)  Elastic  Collisions: 

The  number  of  electrons  per  unit  volume  with  velocities  in  the  range 
dc  ) which  are  scattered  through  an  angle  ( y , & ) into  the  solid 
angle  d(j  by  elastic  collisions  with  molecules  in  the  range  (C,  C + dC 
in  a time  interval  dt  is. 


dn  = f(c)dc  F(C)  dC  g a (y,g)  du  dt 


(6) 


where  g is  the  magnitude  of  the  relative  velocity  and  (y,g)  is  the  differ- 
ential elastic  cross  section.  The  scattering  is  defined  by  a polar  angle 
V and  azimuth  angle  l/>  which  give  the  rotation  of  the  relative  velocity  in 
the  center  of  mass  coordinate  system.  The  functions  f(c)  and  F(C)  are 
the  velocity  distributions  of  the  electrons  and  molecules  respectively.  (The 
spatial  dependence  has  been  suppressed  here  since  all  collisions  are  considered 
point  encounters  involving  no  change  in  position.  ) The  rate  of  change  of 
the  distribution  function  is  found  by  dividing  dn  by  dc  and  dt. 
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Consider  now  a background  gas  which  is  perfectly  cold,  i.  e.  F(C)  = 

Ni(C).  This  is  an  excellent  approximation  in  describing  collisions 
between  electrons  and  heavy  particles  when  the  electron  mean  energy  is 
much  greater  than  the  gas  thermal  energy.  The  relative  velocity,  g,  is 
now  equal  to  the  initial  electron  velocity,  c.  After  integrating  Equation  (6) 
over  C we  have, 
dn 


dc  dt 


f ( c ) Nc<Te(  y,  c)  d<j 


(7) 


where 


Au  = siny  dy  d^ 


Each  collision  term  can  be  divided  into  two  parts,  one  being  the  rate  at 
which  electrons  are  scattered  out  of  a given  element  in  velocity  space  and 
the  other  being  the  rate  at  which  they  are  scattered  in.  The  rate  of  scatter- 
ing out  can  be  found  from  Equation  (7)  by  integrating  over  all  possible 
scattering  angles.  This  gives. 


m 


1(c)  NcQe(c) 


(8) 


e,  out 


where  Q ( c ) = fa  (Y,c)du 

e J e 


The  number  of  electrons  scattered  into  the  interval  (cT*  c +dc)  is  just 
the  sum  of  electrons  scattered  out  of  all  other  intervals  (<T^,  (Tj  + d^J ) 
which  have  final  velocities  in  (c^  c*+dc).  That  is, 


= J’d'cTj  f(<Tj  ) NclJdu>  ae  (y,  c t ) i (c^-cl 


(9) 


e,  in 


where  c^  is  the  final  velocity  of  an  electron  with  initial  velocity  c^.  The 
delta  function  in  Equation  (9)  is  given  in  cylindrical  coordinated  by 


6 (T2-  cl 


6 (c2-c)  6{d?-  0)6  (<t>z-<t>) 


(10) 


c2  sin  ez 
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For  an  elastic  encounter  with  a stationary  target. 


C2  = C 


.[ 


1 - 


2 Mm 
(M+m) 


( i - cosy) J 


1/2 


(ID 


and  c^cos  © ^ 


m 


M 


c.cos  © . + 77 c.  (cos©  cosy  + sin©  siny  cos 

M+m  1 1 M+m  1 1 1 

(12) 


where  m and  M are  the  masses  of  the  electron  and  the  molecule  respectively. 

We  now  change  the  variable  of  integration  in  Equation  ( 9 ) from  c^  to  c 2 
and  expand  all  functions  of  Cj  in  a Taylor  series  around  c^,  since  ci*c2 
is  of  order  m/M  which  is  small.  Thus, 

dc 


:i ' [* +ir  (1  - cosy)  ] dc2 


(13) 


and  CjVi,  0j)  C&(  y,  cp  - c23f(c2>  ©pajy,  c^) 

+ (C1'C2)  ~dZ  [c3f(c’  9\)  ffe(  y,c)]  (14) 

c = c? 


. . m , , ... 

with  c.  - c~  — c0  —(1-cos  / ) 

1 l l M 

Using  these  expressions  in  Equation  ( 9 ) and  integrating  over  c.,  and  4> 
we  have 

IQ  f \ r r i (0  ) 

(Tr)~  N/sin©,d©i  f(c,  ©^c  f du  o (7.  c)  — 

0.  in  2 


(15) 


^ ~§^  c4  J sin  ©jd©L  f(  c,  ©p  Jdcj  (1  -cosY)  ajY.c) 


<J  (©2  -©  ) 

sin  © 


(16) 


The  second  term  on  the  right  is  a result  of  the  small  amount  of  energy  trans- 
ferred from  the  electron  to  the  target  in  the  form  of  recoil  motion, 
b)  Inelastic  Collisions: 

Since  the  rate  of  scattering  out  is  independent  of  the  final  electron  velocity, 
we  have,  in  direct  analogy  with  the  elastic  case. 


(-§7)  = Nf(c.©)  cQ.(c) 

' dt/h,  out  h 


(17 


14 


— 


where  Q (c)  = / O ( y,  c)  d(J  for  each  inelastic  process.  (There  is  an 

h J h 

implied  summation  over  h when  more  than  one  inelastic  process  exists.  ) 


In  an  inelastic  collision  the  final  velocity  is  related  to  the  initial  velocity 
by  , 


1 2 

TmC2 


1 


2 mCl 


(18) 


where  is  the  kinetic  energy  transferred  to  internal  energy  of  the  mole- 
cule. (Here  recoil  has  been  neglected  for  inelastic  collisions.  It  can  be 
included  along  with  the  recoil  term  for  elastic  scattering  if  necessary.) 

From  this  relation  it  follows  that  Cjdcj  = c^dc^.  We  can  now  replace  o&  with 


O-,  in  Equation  (9  ) and  integrate  over  c?  and  <t>  to  obtain, 


fff)  = N /sin^de  1f(C+i  Jdu  ^ (y.c4)_l.(  6 1 

' m,  in  sin  0 


6 ( 6o-0) 
2 


(19) 


The  angles  before  and  after  the  collision  are  related  by, 


c cos  6,  = c^cosOj  cosy  + sinOj  siny  cos  lA  ). 


(20) 


where  y mc^  = y me  + ^ 


c)  Superelastic  Collisions: 


The  superelastic  collision  term  is  obtained  from  the  corresponding  inelastic 
collision  term  by  changing  ^ to  - to  a and  N to  N*.  The 


reverse  cross  section,  <Tg,  is  related  to  the  forward  by  detailed  balance, 


l.  e. 


V y’  Cl> 


<My,  c 


(2i : 


, 1 2 

where  y me ^ 


1 2 
— me 


1 


(22) 
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Using  these  expressions  in  Equations  ( 17 ) and  (19)  we  have. 


i 

A 


N 


s,  out 


*f(c,0)  c+Qh(c+) 


and 


= N*  c ^sin  fljdOj  £(c_  6 J du  o^(Y,c) 


a , in 


6 ( 02  - 0 ) 

sin  $ 


(23) 


(24) 


where 


T mc 


1 2 a 

T“  1 eh 


d)  Ionizing  Collisions: 

First,  we  define  the  probability  per  unit  time  that  an  electron  with  velocity 
c^  making  an  ionizing  collision  produces  an  electron  (scattered  or  ejected) 
with  velocity  in  the  range  (c^,,  c^+dc^)  as, 

NclQi(cl*c2)dc2 

Then,  by  conservation  of  energy,  another  electron  is  also  produced  with 
velocity  c^  given  by, 


1 2 
2 mC3 


Tmci 


1 2 

— me  - €. 

2 2 l 


(25) 


where  is  the  ionization  energy.  (We  have  assumed  for  the  present  that 
the  scattered  and  ejected  electrons  are  produced  isotropically.)  The  collision 
term  for  scattering  out  by  ionizing  collisions  is  then 

c_/y/Z~ 

(ft).  = Nf(c>0)  c f Qi(c’c2)dc2  (26) 

i,  out  JQ 

1 , 12 

where  as  before,  — mc‘  = —me  ± The  upper  limit  of  integration 

“ * L 1 

is  the  velocity  at  which  the  scattered  and  ejected  electrons  have  equal  energy. 

Above  this  value,  we  begin  counting  the  same  events  over  again.. 
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The  collision  term  for  electrons  scattered  into  the  range  (c,  c+dc)  is,  by 
analogy  with  the  inelastic  term, 

00  C1 

(tt)  = “f/sin0!  d$!  f f^,  0^  ~~  Q.tCj.OdCj  (27) 

' ^i.in  J 

r* 


The  integration  is  carried  out  over  all  velocities  which  can  contribute  to 
the  population  at  c. 

e)  Attachment  and  Recombination: 

These  processes  result  in  the  loss  of  an  electron,  so  the  only  effect  on  the 
distribution  function  is, 

= Nf(c,fl)c[Qa(c)  + Qr(c)]  (28) 

a,  r,  out 

Now  the  Boltzmann  collision  term  is  equal  to  the  sum  of  all  contributions 
scattered  in,  minus  the  sum  of  those  scattered  out. 


All  of  the  terms  representing  electrons  scattered  into  the  interval  (c^  c + dc) 
involve  an  integration  over  the  scattering  angle,  Y , the  angle  before  collision, 
6 , and  a delta  function  specifying  the  orientation  after  the  collision  6 . 

JL  4 

For  example,  the  elastic  term  neglecting  recoil  is. 


i(ft)  =nfBin$1deid4lt(c,Ol)fd0  ^(y,  c) 


6[0z-0)6(<rz-sP) 

sin  Q 


(29) 


e,  in 


The  delta  function  can  be  eliminated  by  noting  that  d<Oy=  sin#.,  dO  ^ d <t>^ 
and  performing  the  integration  over  6^  and  0^.  Thus, 

I/2i)  = N/’sind.  d$.  d0,  f(c,  9.)  a (y.c)  (30) 

c \ at  / . J 1 1 1 1 e 

e,  in 


where  cosy 


cos  6 cosQ  + 
1 


sinflj  sin#  cos(0^ -0). 


A similar  transfor- 
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mation  can  be  made  on  the  remaining  scattered- in  collision  terms.  Given 
a differential  cross  section,  ae,  the  integral  in  Equation  (30)  can  be  thought 
of  as  an  operator,  Ke>  acting  on  f.  Then,  if  we  define  an  operator,  K,  as 
the  sum  of  operators  for  each  collision  process,  we  can  write  Equation  (4) 
as 

=7e7SkI  (f)’/2  lK-Q'£,)f 

where  Q(<)  is  the  total  cross  section  for  electrons  with  kinetic  energy, 

(=  t+ri. 


The  operator,  K,  acts  in  the  two-dimensional  space  (£,>?  ),  i.  e.  the  value 
of  Kf  on  the  plane  O = 0 is  independent  of  the  values  of  f for  0 t <i Thus 
K can  be  represented  numerically  as  a four-dimensional  matrix.  Since  the 
range  in  each  dimension  must  be  subdivided  into  at  least  100  intervals,  the 
storage  of  such  a matrix  is  completely  out  of  the  question.  The  functions  Kf 
and  f,  defined  in  terms  of  three  independent  variables,  would  require  three- 
dimensional  arrays  for  their  storage,  also  highly  impractical.  However, 
the  function  f does  not  have  to  be  saved  from  one  iteration  to  the  next  in 
solving  Equation  (31).  Only  Kf  must  be  stored,  and  this  can  be  done  by  making 
a simplifying  assumption.  In  most  practical  cases,  the  function  Kf  is  nearly 
isotropic  in  velocity  space.  Therefore,  we  will  expand  the  angular  depen- 
dence of  Kf  in  Legendre  polynomials,  e.g. 


K f = 
e 


= £ A/(c)  Pf  (cos  6) 
/ = 0 


(32) 


whe  re 


A(  (c) 


2 /+  1 
2 


/*  £(fi)  p'  ,c°s*1  d,cos<>) 

-1  e, in 


(33) 


It  should  be  pointed  out  at  the  beginning  that  the  expansion  in  Equation  (32) 
is  quite  different  from  the  usual  expansion  of  f.  In  particular,  if  0^(V*  c) 
can  be  represented  as  the  sum  of  n Legendre  polynomials,  (cos'/),  then 
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the  sum  in  Equation  (32)  will  also  have  n nonzero  terms,  regardless  of  the 
angular  dependence  of  f.  As  a special  case,  if  o^fc)  is  isotropic,  then  Kf 
is  always  independent  of  0 and  can  be  represented  by  one  term. 

Now,  using  Equation  (29)  in  (33),  we  have 

A^(c)  = Nc  1 J sin0J  d^  f(c,  0j)  J duj  a 2(y,  c)  (cos02)  | 

The  addition  theorem  for  spherical  harmonics  gives 

l 

P/  (COS02)  £ Pl  (cos6l)Y7  (y>)  I 

' m = -/ 

which,  after  integrating  over  the  angle  4> , becomes 
Zjt 

J P(  (cos  #2)  d |//  = 27rp^  (cos  0^ ) P^  (cos  Y)  ( 

0 

When  this  is  inserted  into  Equation  (34),  the  integrals  can  be  separated 
and  Equation  (32)  becomes, 

(■ff)  = Nc  f)  777T  fi(c)  Q/  <c)  pi  <cos^)  (: 

e,  in  „ „ 


whe  re 

Qf(c)  = 4nof(c) 

= 47r-  2 J d(coSy  ) <r(y,  C)  p^  (cosy) 

and 

, , . 2/  + 1 ( 
{1{C)  = 2 7 

d(cos  # ) f(c,  0 ) P(  (cos  0 ) 

In  evaluating  the  integrals  for  the  recoil  term,  there  is  an  additional  fac- 
tor of  (1-cosy)  to  be  considered.  The  integration  over  scattering  angle, 


d(cosy)  (1-cosy ) a (y,  c)  p^  (cosy) 


can  be  carried  out  by  employing  the  recursion  relation, 


to  obtain 


xP 


l (*)  = 


l + 1 
2/  + 1 


P/  +1(X) 


X 


2/  + 1 


P/-  1<*> 


(38) 


2/  + 1 


[at 


(c) 


i±  1 

2^+3 


% 


.(c)  - 


2/  - 1 - 1 


(c) 


The  complete  Boltzmann  collision  term  can  now  be  written, 

(ff)  ‘ Nc£  TiTT 

c 1=0 


m N a 4 V'  1 r / x IX-,  x 
M c2  ac  c X Zl  + 1 [f(c)  [°£(c)  - 


2/ 
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It  has  been  assumed  that  secondary  electrons  produced  in  ionizing  colli- 
sions are  emitted  isotropically.  The  contributions  from  superelastic 
collisions  have  been  omitted.  Because  of  the  very  small  ratio  of  electron 
mass  to  neutral  mass,  the  recoil  term  can  also  be  neglected  for  most 
electric  discharge  applications. 


The  Boltzmann  equation,  in  terms  of  the  electron  energy,  € = — , is  then 
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The  integro-differential  equation  above  can  be  solved  for  the  distribution 
function  f (£,?}, <p),  or  given  a nonuniform  E/N(<p)  and  a known  set 

of  cross  sections.  The  program  which  was  developed  to  do  this  will  be 
described  in  the  following  section. 


4.  2 Electric  Field  and  Ion  Kinetics 

Poisson's  equation  for  the  scalar  potential,  defined  by  E = i£ 


= 77  K - n 


ne) 
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where  n+,  n , and  n^  are  the  positive  ion,  negative  ion,  and  electron  num- 
ber densities,  respectively.  For  the  present,  only  one  type  of  positive 
ion  and  negative  ion  will  be  considered.  The  equations  can  easily  be  modi- 
fied to  include  more  than  one. 


The  heavy  ions  are  governed  by  the  moment  equations,  which  represent 
conservation  of  number,  momentum  and  energy.  The  continuity  equations 
are 
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where  v and  v are  the  positive  and  negative  ion  drift  velocities.  The 

terms  on  the  right  hand  side  represent  production  and  loss  due  to  external 

ionization,  S,  Townsend  ionization,  z,  attachment,  a,  detachment,  d,  and 

ion-ion  and  electron-ion  recombination,  r.  and  r , respectively.  The 

l e 

coefficients  involving  electron  impact  are  defined  in  terms  of  the  electron 
energy  distribution  and  the  appropriate  cross  sections,  i.e. 
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Notice  that  these  coefficients  are  functions  of  position,  as  well  as  of  the  bound- 


ary conditions  which  go  into  calculating  f^ff,'/7). 


The  mean  velocities,  v^  and  v , are  determined  by  momentum  conservation, 


n (v  • V)  v + V * p - en  E = -m  n v (v  - V)  - m v £ 
- -r+  ,-f  r+  + + + -V  + 


(44) 


m n (v  • V)  v + V • p + en  E = -m  n v (v  - V)  - m v S 


(45) 


where  m and  m are  the  masses  of  positive  and  negative  ions,  respectively 


and  V and  V are  the  total  collision  frequencies  of  these  ions  with  neutrals. 

+ 


The  mean  neutral  velocity  is  V.  An  equation  of  state  relates  the  ion  pressure, 


p^,  to  the  mean  thermal  velocity,  c^,  and  the  temperature,  T±  . Thus 


1 / 3m+n  tc± 


n ± kT± 


(46) 


The  series  of  moment  equations  for  heavy  ions  is  terminated  by  assuming 
their  thermal  motion  to  be  in  equilibrium  at  the  local  gas  temperature,  T. 
This  approximation  is  fairly  good  because  of  the  nearly  equal  masses  of 
ions  and  neutrals  and  the  fact  that  at  high  pressures  the  sheath  is  collision 
dominated. 


Let  us  consider  now  a plasma  consisting  of  positive  ions  and  electrons  in 
a one-dimensional  geometry.  Then  using  Equation  (42)  in  (44)  without  gas 
flow, 
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If  diffusion  is  negligible  compared  to  inertia  (n  — const)  and  the  drift  energy 

2 + 

is  large  compared  to  thermal  energy  (m  v » kT  ),  then 
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or  in  terms  of  the  potential,  (fi, 


d t P 


2E 


2m . 


v. 


(49) 
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The  collision  frequency  for  momentum  transfer  is  given  by 
V.  = NQv, 


(50) 


where  Q is  dominated  by  charge  exchange  and  can  be  considered  independent 
of  velocity. ^ 


Given  E(</>),  Equation  (29)  is  integrated  to  obtain  v+(<P).  Then  the  positive 
ion  number  density  is  determined  from 


ev 


J-h 


ev 


(51) 


where  j is  given  and  is  obtained  from  the  solution  of  Boltzmann's 

equation  for  electrons.  Finally,  Equation  (41)  is  written  in  terms  of 

-k  (-§-)  = t hw,)  • (52) 

and  integrated  to  yield  a new  E^\c/)).  This  is  used  in  Equation  (29)  and 
the  procedure  is  repeated  until  convergence  is  reached.  The  Boltzmann 
equation  for  electrons  must  be  solved  for  each  new  field,  E^. 


5.0 


COMPUTER  PROGRAMS 


Two  computer  programs  are  described  briefly  in  this  section.  The  first  was 
developed  entirely  under  the  present  contract  to  solve  the  two-dimensional 
Boltzmann  equation  for  electrons  which  is  then  coupled  self-consistently  with 
the  ion  kinetics,  gas  dynamics,  and  field  equations.  The  second  is  a program 
for  solving  the  boundary  layer  equations  along  a flat  plate  which  was  modified 
in  a previous  contract  to  include  a heat  addition  term  defined  by  the  electric 
discharge. 

5.1  Program  SHEATH 

The  approach  used  to  solve  the  Boltzmann  equation  involves  a numerical 
integration  along  the  characteristics,  £ - C = const,  combined  with  an  itera- 
tion scheme  for  evaluating  the  collision  integral,  Kf.  Initially  it  is  assumed 
that  there  is  no  contribution  from  scattered-in  electrons,  i,  e.  Kf  = 0,  and 
Equation  (40)  is  integrated  in  both  the  forward  and  backward  directions  for 
each  value  of  rj.  This  gives  a first  approximation  to  the  distribution  function 
f^  ^ from  which  a new  Kf^  is  calculated,  and  the  procedure  is  repeated 
until  convergence  is  reached. 

The  operations  required  to  generate  successive  approximations  to  the  col- 
lision integral  are  performed  as  the  integration  proceeds,  so  that  the 
function  f does  not  have  to  be  stored  between  iterations.  Only  Kf  has  to 
be  retained.  By  performing  the  integrations  in  a specific  order,  the  same 
array  can  be  used  to  store  the  old  and  the  new  values  of  Kf.  If  Kf  is  repre- 
sented by  the  first  two  terms  in  a Legendre  expansion,  then  two  two-dimen- 
sional arrays  are  needed.  These  arrays  are  trapezoidal  in  shape,  since 

f(^,r7,i.')  = 0 for  £ ><P+$  , where  £ is  the  maximum  value  of  longitudinal 

o o 

energy  for  electrons  leaving  the  cathode. 

The  storage  arrays  for  Kf  are  arranged  in  core  as  in  Figure  1 to  minimize 
central  memory  requirements.  The  dimensions  were  chosen  to  give  1 eV 
resolution  for  a 200  V cathode  fall.  The  energy  distribution  of  electrons 
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emitted  from  a cathode  by  positive  ion  impact  has  a cutoff  at  f - 2t 0 , where 

i w 

*s  the  ionization  potential  of  the  gas  and  y?w  is  the  work  function  of  the 
cathode  material.  This  cutoff  is  usually  less  than  20  eV,  so  that  20  mesh 
points  at  the  cathode  are  again  sufficient  for  1 eV  resolution.  The  area  of 
the  array  in  Figure  1 is  48,  240jq  (136,  160g)  words.  Obviously,  the  storage 
of  f (or  of  Kf  without  angular  expansion),  which  would  require  at  least  ten 
times  this  much  core,  is  impractical. 


Boundary  conditions  are  specified  at  the  cathode  and  at  the  anode.  At  the 
cathode,  the  energy  and  angular  distribution  of  emitted  electrons  are  given 
At  the  anode,  the  distribution  of  backscattered  and  ejected  secondary  elec- 
trons is  given  in  terms  of  the  incident  primary  flux.  At  the  turning  points, 
the  forward  and  backward  fluxes  are  equal. 


The  algorithm  used  for  the  numerical  integration  of  Equation  (40)  is  the 
following. 

■ [cuitCi  + (i7-^77)  <«> 

. 1 k vl/2 

where  C - eE/N  \J~ ) and  ^ is  the  cncrgy  step.  The  scheme  is  second 

order  and  absolutely  stable.  The  relative  error  at  each  step  is, 

(54) 
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The  error  is  small  as  long  as 

< 0.  5 E/NQ 

This  condition  is  easily  satisfied  near  the  cathode  surface,  where  E/NQ 
and  fj  is  typically  20  eV.  it  can  also  be  satisfied  in  the  positive 
column,  where  E/NQ  - 1 eV.  However,  in  the  negative  glow,  the  field  be- 
comes very  small  and  may  even  reverse  directions.  Under  these  conditions, 
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the  algorithm  of  Equation  (53)  cannot  be  used  and  in  fact  the  choice  of  vari- 
ables in  Equation  (40)  is  no  longer  appropriate,  since  position  becomes  a 
multivalued  function  of  the  potential. 

In  the  region  where  E/NQ«1  eV,  the  Boltzmann  equation  can  better  be 
expressed  as, 


1/2 

Ji  + ®E(y)  f|  = N(y)  ( f ) )k-Q(£)|  < 


(55) 


Now  a fixed  position  step,  4y,  is  chosen  with  the  condition  that, 
dy<  0.  283/NQ,  where  1/NQ  is  the  electron  mean  free  path.  This  choice 
of  variables  cannot  be  used  in  the  high  field  region  near  the  cathode.  There- 
fore, to  fully  define  f in  every  region  of  the  discharge,  we  must  start  with 
coordinates  ($,r?,V>)  in  the  cathode  fall,  change  to  coordinates  (^,  »7,  y)  m 
the  negative  glow  and  Faraday  dark  space,  and  back  to  (?- rl>^  ) m the  P°s1' 
tive  column.  Basically  the  choice  of  coordinates  depends  on  whether  the 
distribution  function  is  changing  more  rapidly  with  respect  to  position  or 

potential. 

With  the  present  program,  only  the  cathode  fall  region  has  been  examined. 
The  code  is  currently  set  up  for  integration  in  the  plane  V = 0,  i.e.  electrons 
are  restricted  to  moving  either  with  or  directly  against  the  electric  field. 
This  simplification  reduces  the  running  time  by  about  a factor  of  ten  and 
allows  rapid  checkout  and  modification  of  the  program.  The  distribution 
function  is  normalized  in  a one-dimensional  velocity  space,  i.e. 


oo 

^ff+(£,^)+f  )]‘  ~1/2d(  = ne(C>) 


(56) 


where  f is  the  distribution  of  electrons  moving  toward  the  anode  and  f_  is 

the  distribution  moving  toward  the  cathode.  These  functions  satisfy  Equation(31 

with  one  of  the  following  pairs  of  collision  terms. 
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(a)  Forward  scattering: 
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(b)  Isotropic  scattering: 
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(c)  Anisotropic  scattering 
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The  cross  sections  for  attachment  and  recombination  are  small  at  high 
energies  and  are  therefore  not  included  in  the  analysis  of  the  cathode  fall 
region.  They  will  be  important,  however,  in  the  negative  glow,  Faraday 
dark  space,  and  positive  column.  The  secondaries  produced  by  ionizing 
collisions  are  assumed  to  be  evenly  distributed  in  energy  from  zero  to 
f-  e ■ . This  is  a good  approximation  near  threshold.  However,  at  higher 
energies  the  distribution  is  peaked  near  zero  and  ( - t ■,  with  a minimum 
at  (<  - ( J/2.  Previous  workers  have  found  that  the  actual  distribution  of 
secondaries  makes  little  difference  to  the  overall  results.  This  aspect  may 
be  investigated  with  the  present  code  in  the  future. 

The  input  parameters  required  by  the  program  are  (1)  the  electric  field 
strength  at  the  cathode,  (2)  the  reduced  gas  pressure  (this  is  the  actual 
pressure  reduced  to  standard  temperature,  i.e.  p^  = p(273.  1 6°K /T ),  and 
(3)  the  voltage  difference  between  cathode  and  anode.  The  function  E/N(<^) 
is  then  generated  internally  for  either  a constant  or  linear  dependence  on  y. 
The  cross  sections  are  read  in  next,  following  a standardized  format  used 
previously  by  the  author  in  an  equilibrium  Boltzmann  code.  (The  specific 
input  deck  structure  will  be  provided  in  the  user's  manual.  ) And  finally, 
the  energy  distribution  of  electrons  leaving  the  cathode  is  given. 

The  output  data  supplied  by  the  program  may  be  varied  to  suit  the  users' 
needs.  Presently,  the  listing  is  set  up  as  follows:  The  input  parameters 
and  cross  sections  are  printed  first.  Then,  in  order  to  monitor  the  con- 
vergence, the  electron  current,  je,  as  a function  of  potential, and  the  energy 
distribution,  f^,  at  the  anode  are  listed  after  each  iteration.  The  criterion 
for  convergence  is  that  the  electron  current  at  the  anode  change  by  no  more 
than  one  percent.  When  this  is  reached,  the  following  quantities  are  printed 
as  functions  of  (fi:  the  electron  mean  energy,  e . the  distance  from  the 

cathode,  y,  the  electron  number  density,  ne,  the  electron  drift  velocity, 

1 6 j 0 

v , Townsend's  first  ionization  coefficient,  defined  as  a - -J- — , and 

e Je  dY 

E/N.  Even  though  the  entire  distribution  function  f (^ , 17  , c/? ) is  never  known 
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at  any  given  time  in  the  calculation,  it  can  be  printed  out  if  desired.  All 
the  important  information,  however,  can  be  extracted  by  listing  f(£,t7,</?0) 
at  specific  locations,  tpQ,  and  by  listing  various  moment  integrals,  such  as 
drift  velocity,  mean  energy,  diffusion  coefficient,  etc.,  as  functions  of 
position. 

5.  Z Program  CATFAL 

CATFAL  is  a computer  program  for  calculating  laminar  and  turbulent  boundary 
layer  development  in  compressible  flow  including  the  effects  of  cathode  fall. 

The  program  is  an  extension  of  a program  written  at  Princeton  University 
and  was  developed  under  contract  with  Air  Force  Weapons  Laboratory  to 
analyze  a high  power  CO  CW  gas  laser.  The  laser  configuration  is  illustrated 
in  Figure  2.  Gas  flows  through  a two-dimensional  supersonic  nozzle  at  Mach 
3.  5 to  4.  0.  A high  voltage  (2-5  kV)  is  applied  perpendicular  to  the  flow  and  a 
discharge  established  in  the  region  where  the  optical  power  is  to  be  extracted. 

The  electrical  discharge  is  stabilized  by  means  of  an  electron  beam  whose 
electrons  have  an  energy  of  approximately  150  keV.  In  such  a laser,  a large 
potential  gradient  exists  in  the  proximity  of  the  cathode  which  is  usually  refer- 
red to  as  the  "cathode  fall  region". 

In  this  region,  due  to  the  large  potential  drop,  strong  electric  fields  are 
established.  This  results  in  a substantial  amount  of  power  in  the  form  of 
heat  being  dissipated  in  the  boundary  layer  next  to  the  cathode.  Such  action 
produces  a disturbance  which  propagates  into  the  flow  region.  Because  of 
the  stringent  requirements  of  flow  homogeneity  in  CO  CW  lasers,  the  CATFAL 
computer  program  was  written  to  estimate  the  magnitude  of  the  flow  dis- 
turbance. The  code  predicts  the  effect  of  cathode  fall  on  the  growth  rate  of 
the  boundary  layer.  The  analytical  model  for  CATFAL  was  developed  by 
Claudio  Parazzoli,  formerly  of  Northrop.  Experiments  were  performed  to 
verify  the  numerical  predictions  of  CATFAL  and  those  results  are  published 
in  Reference  1 . 
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5.  2.  1 Analytical  Model 


In  this  section  we  shall  briefly  describe  the  boundary  layer  equations  for  a 
turbulent  compressible  flow  on  a flat  plate  and  we  shall  put  them  in  a form 
manageable  for  numerical  calculations. 


Boundary  Layer  Equations 

The  equations  governing  the  flow  of  a compressible,  two-dimensional  boundary 
layer  are: 


Continuity  Equation 
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Momentum  Equation 
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Energy  Equation 
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Figure  3 illustrates  the  relevant  quantities  in  Equations  (60)  through  (62). 


Equations  (60)  to  (62)  apply  both  to  laminar  and  turbulent  boundary  layer 
provided  that  the  definitions  for  r and  q are  taken  to  be 
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where  -u'  v'  is  the  kinematic  Reynold's  stress;  -v'  h‘  is  the  kinematic  turbulent 

U 

heat  flux;  u = is  the  molecular  kinematic  viscosity,  and  v = — is  the 


g PCr 


molecular  kinematic  conductivity. 


Then  we  define  an  effective  viscosity  and  an  effective  conductivity  as: 
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Finally,  the  following  quantities  are  introduced  for  convenience  in  the  calcula- 
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where  <5  is  the  boundary  layer  displacement  thickness  given  as: 

CO 

6*  = f (Pe  U(x)  - pu  (x,  y) ) / Pe  U(x)  dy  (67) 

o 

When  written  in  terms  of  the  new  variables  defined  in  Equations  (64)  to  (66), 
Equations  \ 60 ) to  (62)  become: 
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where  1(d)  is: 


J >(x)  - 


with  x = - p and  Vis  the  variable  specific  heat  ratio  across  the  boundary 
layer.  Equations  (68)  to  (70)  are  similar  to  expressions  II - 1 4 , 11-15 
and  11-16  in  Reference  1.  The  major  differences  are  in  Equations  (69)  and 
(70).  Equation  (69)  contains  the  heat  addition  term  in  the  R.H.S. 


which  is  absent  in  Equation  11-15  of  Reference  1.  Equation 

P U (h°  - h ) 
e e r 

(70)  accounts  for  a variable  specific  heat  ratio  across  the  boundary  layer,  and 
reduces  to  Equation  11-16  of  Reference  1 if  we  take  f - "Ye  . In  the  derivation 
of  Equation  (70),  the  specific  heat  ratio  of  ye  , in  the  free  stream,  has  been 
taken  to  be  constant. 

References  1 and  47give  detailed  discussions  of  the  form  of  normalized  effective 


vi  scosity  *»e  - 


— , and  effective  heat  conductivity  feg  - 


5,  z.  2 Computer  Code 

The  original  gas  dynamics  program  was  written  by  H.  J.  Herring  and 

47 

G.  L.  Mellor  at  Princeton  University.  The  program  performs  a numerical 
integration  of  the  equations  of  motion  for  a compressible  two-dimensional 
boundary  layer.  Calculations  may  be  carried  out  for  both  laminar  and 
turbulent  flow  for  arbitrary  Reynolds  number  and  freestream  Mach  number 
distribution.  Planar  or  axisymmetric  bodies  with  wall  heating  or  cooling, 
wall  suction  or  blowing  and  rough  or  smooth  wall  can  be  included  . 

Modifications  were  made  to  the  program  to  include  bulk  heat  addition  due 
to  the  cathcde  fall  of  a glow  discharge  in  the  two-dimensional  planar  case. 
Boundary  layer  equations  were  modified  to  include  a heat  addition  term. 
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The  program  can  handle  arbitrary  gas  mixtures  of  up  to  four  components. 
Calculation  of  molecular  viscosity  and  heat  conductivity  of  gas  mixtures  as 
a function  of  temperature  is  included.  Allowance  of  variable  specific  heats 
and  ion  and  electron  drift  velocities  in  a gas  mixture  are  computed  as  a 
function  of  E/N  and  tabulated  for  utilization  in  computing  the  power  term. 
Poisson's  equation  and  current  continuity  equations  are  then  solved  using 
iterative  techniques  and  Runge-Kutta  integration  to  satisfy  certain  boundary 
conditions. 


The  MAIN  program  was  modified  to  compute  viscosity  and  heat  conductivity 
for  an  arbitrary  gas  mixture  of  up  to  four  components.  Subroutine  VISC 
provides  fourth  order  polynomial  curve  fits  of  viscosity,  heat  conductivity 
and  Prandl  number  as  a function  of  temperature.  SHRSUB  then  computes  the 
specific  heat  ratio  for  the  gas  mixture. 

The  MAIN  program  iterates  to  find  velocity  (FP)  and  enthalpy  (GP)  profiles 
within  desired  convergence  criteria.  Before  each  FP-GP  ite ration,  sub- 
routine POWER  solves  Poisson's  equation  and  the  current  continuity  equa- 
tions using  a fourth  order  Runge-Kutta  approximation.  An  iterative  procedure 
is  used  to  satisfy  the  boundary  conditions.  Power  dissipation  :n  the  vicinity 
of  the  cathode  is  computed  and  included  in  the  boundary  layer  calculation. 

The  program,  written  entirely  in  FORTRAN,  has  been  run  on  the  CDC  175  and 
has  a core  requirement  less  than  160,000  words  (octal).  Typical  running 
time  with  heat  addition  is  normally  3-4  minutes.  Output  at  each  X position 
is  written  on  file  9 and  can  be  saved  on  tape  or  permanent  file  for  subsequent 
restart  if  desired. 
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NUMERICAL  RESULTS 


6.  0 

The  computer  program,  SHEATH,  was  developed  in  stages  over  a period  of 
time  and  was  checked  at  each  stage  by  making  comparisons  with  eases  pub- 
lished in  the  literature.  All  of  these  runs  were  made  assuming  constant  gas 
density  and  a known  electric  field  variation.  One  self-consistent  field  solution 
is  given  at  the  end  of  Section  6.  1.  The  difficulties  encountered  in  extending 
the  self-consistent  solution  into  the  negative  glow  are  discussed  in  Section  6.2. 


6.  1 Nonequilibrium  Electron  Kinetics 


The  first  runs  were  made  in  mercury  for  comparison  with  the  Monte  Carlo 
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simulations  of  Parker,  et  al.  The  ioni/.ation  cross  section  of  Kieffer,  et  al 
was  fitted  to  the  functional  form, 
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where  A = 5.  25  C <7max>  - 10.4  eV,  and  £Tnvix  = 6.  5A1".  (At  this  point 
the  subroutine  for  reading  cross  section  data  in  numerical  form  had  not 
been  incorporated  into  the  program.)  Forward  scattering  was  assumed 
with  no  backscatter  from  the  anode.  The  electrons  were  allowed  to  fall 
300  V at  a constant  E/pQ  of  4000  V/cm/Torr.  The  energy  distribution  of 
electron  current  striking  the  anode  is  shown  in  Figure  4 . 


A similar  case  was  examined  in  Reference  15  with  a gap  separation 
of  1000  V.  The  two  distributions  compare  very  well  below  about  250  eV, 
indicating  that  the  low  energy  electrons  have  reached  an  equilibrium  with 
the  field.  A plot  of  f a /f  indicates  that  the  distribution  is  very  nearly 
Maxwellian  in  this  region, with  a temperature  of  ~ 1 00  eV.  The  peak  in  the 
energy  distribution  above  300  eV  represents  those  electrons  which  have  suf- 
fered no  ionizing  collisions.  The  5 eV  width  of  this  peak  is  a reflection  of 
the  energy  distribution  leaving  the  cathode.  The  gap  below  300  V has  a 
width  of  - W,  where  W is  the  width  of  the  primary  peak.  As  pointed  out 
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earlier,  W = <•  - 2</>w>  so  the  width  of  the  gap  is  just  equal  to  twice  the  work 
function  of  the  cathode  material,  ip w> 

The  electron  current  density  and  first  Townsend  ionization  coefficient  are 
shown  in  Figure  5.  There  is  no  ionization  until  the  electrons  have  first 
been  accelerated  to  the  ionization  threshold.  This  occurs  at  y = €^/ E = 

0.0104  cm.  The  Townsend  coefficient  rises  rapidly  above  threshold  and 
overshoots  the  eventual  equilibrium  value  of  aQ  — 4.  0 cm  (This  is  to 
be  compared  with  the  value  obtained  in  Reference  15  of  0o~  5.  1 cm 
at  the  same  E/N  and  pressure,  but  using  a different  ionization  cross  section.  ) 
The  nonequilibrium  behavior  of  a over  most  of  the  discharge  gap  is  obvious 
from  the  figure.  However,  the  multiplication  factor,  M = je(d)/je(o)  = 3.6, 
is  very  close  to  the  value  e a°^  = 3.3  which  would  obtain  for  constant  a =Q0. 

The  foregoing  example  was  meant  to  simulate  conditions  typical  of  dielectric 
breakdown  in  low  pressure  gases,  where  the  electric  field  is  uniform.  In 
order  to  investigate  the  nonequilibrium  behavior  of  electrons  in  the  cathode 
region  of  a glow  discharge,  we  examined  the  following  problem.  Let  the 
electric  field  fall  linearly  from  1000  V/cm  at  the  cathode  to  50  V/cm  at  the 
midplane  and  remain  constant  from  there  to  the  anode.  The  reduced  pres- 
sure is  0.  25  Torr  and  the  potential  difference  is  again  300  V (270  V in  the 
"cathode  fall"  and  30  V in  the  "negative  glow").  The  results  are  shown  in 
Figure  6 . 

The  behavior  of  a over  the  first  0.  4 cm  is  very  similar  to  that  of  Figure  5, 
even  though  in  the  present  case  the  electric  field  has  fallen  to  less  than 
300  V/cm.  This  is  because  most  of  the  ionization  in  this  region  comes  from 
the  primary  group,  whose  energy  is  f =c’,  independent  of  the  local  field. 

The  Townsend  ionization  coefficient  for  a monoenergetic  group  of  electrons 
is  Q (< ) = NQ*(f  )>  which  has  a maximum  value  of  5.  75  cm  at  60  eV.  This 
maximum  is  never  reached  in  practice,  however,  because  secondaries  are 
produced  at  low  energies  with  lower  ionization  efficiencies  than  the  primary 
group. 
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Figure  5 Current  Density  and  Ionization  Coefficient 
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The  persistence  of  ionization  after  the  electric  field  is  reduced  is  another 

nonequilibrium  feature  evident  in  Figure  6.  The  rate  at  which  Q decays 

from  one  equilibrium  value  to  a lower  one  must  clearly  depend  on  the  degree 

of  angular  scattering.  In  the  present  case  of  forward  scattering  (which  is  actually) 

no  scattering  at  all)  the  persistence  is  most  pronounced.  As  we  shall  see 

later,  in  the  case  of  isotropic  scattering  the  persistence  is  reduced.  The 

situation  with  real  gases  in  the  negative  glow  is  closer  to  the  case  of  forward 

scattering,  since  the  differential  cross  sections  for  energies  above  20  eV 

are  very  strongly  peaked  in  the  forward  direction.  In  fact,  the  beam-like 

character  of  electrons  leaving  the  cathode  fall  and  persisting  through  the 
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negative  glow  has  been  well  established  experimentally. 


The  electron  current  and  number  density  are  also  shown  in  Figure  6. 

Notice  that  je  increases  by  about  a factor  of  four  in  the  negative  glow  or 

drift  region.  With  previous  theories  of  the  cathode  fall  using  equilibrium 

coefficients,  current  growth  in  this  region  would  be  negligible.  Thus 

predicted  current  densities  as  a function  of  cathode  fall  potential  would 
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be  low.  This  is  clearly  evident  in  a comparison  of  Ward's  results  with 
experimental  data. 

As  electrons  are  accelerated  away  from  the  cathode  surface,  the  number 

density  drops  initially  in  order  to  maintain  current  continuity.  Then  as 

ionization  picks  up  the  density  goes  through  a minimum  and  begins  to  rise 
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exponentially.  This  same  behavior  was  observed  by  Tran  Ngoc,  et  al 
in  their  Monte  Carlo  simulation.  The  initial  drop  is  less  severe  if  elec- 
trons are  emitted  from  the  cathode  at  higher  energy,  since  they  reach 
ionization  threshold  sooner.  The  electron  density  is  rising  faster  than 
the  current  at  y = 0.  4 cm  because  the  drift  velocity  is  falling  with  the  field. 

The  next  set  of  runs  was  made  in  argon  with  isotropic  scattering  and  one 
inelastic  process  in  addition  to  ionization.  The  results  are  compared  with 
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the  Monte  Carlo  simulations  of  Sakai,  et  al  in  a uniform  field.  The  cross 

sections  used  here  are  shown  in  Figure  7.  They  are  similar,  though  not 

identical  to  those  employed  in  Reference  5.  The  inelastic  cross  section, 

which  is  meant  to  represent  the  sum  of  all  electronic  excitations,  was 
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taken  from  Schaper,  et  al.  For  elastic  scattering  we  used  the  momentum 
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transfer  cross  section  of  Spencer,  et  al  and  for  ionization,  the  measure- 
ments of  Rapp,  et  al.  ^ 

The  drift  velocity  and  ionization  coefficient  for  E/N  = 565  Td  {1  Td  = 

10  ^ V cm^)  are  shown  in  Figure  8 as  a function  of  distance  from  the 
cathode.  The  most  striking  feature  of  these  curves  is  the  series  of  dis- 
continuities in  the  drift  velocity  beginning  near  the  cathode  and  diminish- 
ing in  amplitude  as  y increases.  Similar  oscillations  had  been  observed 
previously  by  Sakai,  et  al1^  and  were  ascribed  to  the  onset  of  inelastic 
thresholds.  This  behavior  can  now  be  more  fully  explained  as  follows. 

Below  the  first  inelastic  threshold  at  11.5  eV  there  are  only  primary  elec- 
trons which  have  suffered  no  inelastic  collisions.  Due  to  elastic  scattering 
however,  these  electrons  are  moving  both  with  and  against  the  electric  field. 
Their  drift  velocity  is  towards  the  anode  and  increases  as  they  are  acceler- 
ated by  the  field. 

Above  the  first  inelastic  threshold,  there  appears  a group  of  electrons  which 

have  suffered  one  inelastic  collision.  They  have  an  energy  e = <P-  . 

Those  electrons  which  are  moving  toward  the  cathode  will  be  decelerated 

by  the  field  until  they  reach  zero  velocity  at  the  point  y = € /E.  Thus, 

h 

immediately  above  the  first  inelastic  threshold,  there  is  a group  of  electrons 
with  zero  drift  velocity  superimposed  on  the  primary  group.  This  gives  rise 
to  the  first  discontinuity.  The  drift  velocity  then  continues  to  increase  as 
both  groups  of  electrons  are  accelerated  toward  the  anode. 
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A similar  drop  in  the  drift  velocity  occurs  at  the  ionization  threshold  of 
16  eV,  except  that  here  the  energy  loss  is  not  discrete  but  is  distributed 
over  a range  of  values.  A series  of  discontinuities  occurs  at  two,  three, 
and  four  times  the  inelastic  and  ionization  thresholds,  caused  by  electrons 
which  have  undergone  two,  three,  and  four  inelastic  or  ionizing  collisions, 
and  so  on.  The  amplitude  of  the  oscillations  diminishes  as  the  number  of 
secondary  electrons  increases.  These  have  a broad  distribution  of  ener- 
gies and  the  effect  is  smeared  out. 

Also  shown  in  Figure  8 is  the  variation  in  drift  velocity  when  forward  scat- 
tering is  assumed.  There  are  no  oscillations  in  this  case  because  all  the 
electrons  are  moving  toward  the  anode.  The  equilibrium  value  of  drift 
velocity,  which  is  reached  at  about  y = 0.  2 cm,  is  an  order  of  magnitude 
larger  than  for  isotropic  scattering.  This  implies  an  inverse  relationship 
for  the  electron  densities  in  the  two  cases,  since  the  current  densities  are 
almost  the  same  when  normalized  to  the  current  density  at  the  cathode. 

An  important  point  should  be  made  at  this  time  with  regard  to  the  emission 
of  electrons  from  the  cathode.  Experiments  which  measure  secondary 
electron  ejection  by  ion  impact  are  usually  carried  out  in  a high  vacuum 
environment.  In  the  presence  of  a gas  at  high  pressure,  the  yield  can  be 
considerably  reduced.  The  electron  current  density  at  the  cathode  is  not 
only  a function  of  the  rate  of  ejection  but  also  depends  on  the  surface  E/N 
and  on  the  degree  of  scattering  by  the  gas.  In  the  cases  considered  above, 
electrons  were  ejected  from  the  cathode  at  the  rate  of  y per  unit  area.  The 
current  density  in  the  case  of  forward  scattering  was  e7,  while  for  isotropic 
scattering  it  was  only  0.  3 eY  because  70%  of  the  ejected  electrons  were 
backscattered  in  the  gas  and  reabsorbed  by  the  cathode. 

The  Townsend  a coefficient  in  Figure  8 is  zero  below  the  ionization  threshold 
then  rises  rapidly  to  a peak  and  undergoes  a few  oscillations  analagous  to 


48 


those  in  the  drift  velocity.  The  equilibrium  value  of  4.  5 cm  at  y = 0.  2 cm 

12  - 1 

compares  well  with  the  value  obtained  by  Tagashira,  et  al  of  4.3  cm 

7 7 

The  drift  velocity  at  this  point  is  2.  7 x 10  cm/sec  compared  to  3.  2 x 10  cm/eec 
in  Reference  4 . The  increase  in  drift  velocity  from  y = 0.  2 cm  to  0.  5 cm 
is  a consequence  of  the  boundary  condition  at  the  absorbing  anode.  The 
flux  of  electrons  moving  toward  the  cathode  decreases  from  its  equilibrium 
value  to  zero  at  the  anode.  At  this  point  all  electrons  are  moving  in  the 
forward  direction  so  the  drift  velocity  is  the  same  as  that  with  forward 
scattering. 

The  distribution  of  electron  current  striking  the  anode  is  shown  in  Figure  9. 

The  most  prominant  feature  of  this  plot  is  the  series  of  equally  spaced 
peaks  below  100  V.  These  represent  the  groups  of  electrons  which  have 
undergone  0,  1,  2,  . . . inelastic  energies  leaving  the  cathode.  In  this  case, 
a box  of  0.  5 eV  width  has  been  assumed.  The  series  of  plateaus  represent 
electrons  which  have  undergone  one  ionizing  and  0,  1,  2,  . . . inelastic  col- 
lisions. Electrons  experiencing  more  than  one  ionizing  collision  and  all 
secondary  electrons  make  up  the  body  of  the  distribution. 

In  Figure  10  the  electron  energy  distribution  is  plotted  on  a linear  scale  for 
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comparison  with  the  Monte  Carlo  results  of  Sakai,  et  al.  Although  there  is 
considerable  statistical  scatter  in  the  Monte  Carlo  data,  some  areas  of  simi- 
larity and  difference  can  be  pointed  out.  The  Monte  Carlo  results  fall  con- 
sistently below  the  present  results  at  energies  below  10  eV  and  consistently 
above  10  eV  and  20  eV.  The  high  energy  tail  of  the  distribution  is  in  fairly 
good  agreement.  The  differences  may  be  due  to  the  different  cross  sections 
used  or  they  may  reflect  the  approximations  made  in  angular  scattering. 

The  next  run  in  argon  was  made  with  an  electric  field  which  varied  linearly 
with  position,  falling  from  600  V/cm/Torr  at  the  cathode  to  30  V/cm/Torr 
at  the  anode.  Figure  11  shows  the  variation  of  electron  current  and  number 
density  with  distance  from  the  cathode.  The  anode  is  at  0.33  cm  Torr  and 
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Figure  11  Electron  Current  and  Number  Density 
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was  taken  to  be  either  perfectly  absorbing  or  to  reflect  50%  of  the  incident 
electrons.  The  electron  current  is  constant  below  the  ionization  threshold 
which  occurs  at  0.028  cm  Torr.  It  then  rises  smoothly  and  begins  to  level 
off  before  rising  again  in  front  of  the  anode.  The  number  density  behaves 
just  as  the  inverse  of  the  drift  velocity  at  low  energies  since  j = eneW. 

The  falloff  in  approaching  the  anode  indicates  a rising  drift  velocity.  Notice 
that  with  50%  backscatter,  the  effect  of  the  anode  on  the  electron  distribu- 
tion is  reduced. 

The  Townsend  ionization  coefficient  for  this  case  is  shown  in  Figure  12. 

The  form  is  quite  different  from  what  would  be  obtained  by  assuming  an 
equilibrium  energy  distribution  at  the  local  E/N.  The  ionization  in  the 
latter  case  would  be  greatest  at  the  cathode  and  fall  off  as  shown  by  the 
dotted  line.  The  actual  ionization  is  zero  near  the  cathode  due  to  the  thres- 
hold effect  and  is  higher  than  the  equilibrium  value  in  the  low  field  region 
because  of  the  time  required  for  fast  electrons  to  slow  down.  Note  once 
acain  the  effect  of  changing  the  boundary  condition  at  the  anode. 

The  three  functions  which  will  couple  directly  to  the  ion  kinetics  and  the 
field  equations  in  the  self-consistent  model  of  the  cathode  region  arc  the 
ionization  coefficient  and  the  electron  current  and  number  densities  versus 
position.  The  product  of  the  first  two  of  these,  Qje,  gives  the  ionization 
source  density  which  appears  in  the  positive  ion  continuity  equation.  The 
electron  number  density  appears  in  the  Poisson  equation.  The  numerical 
formulation  of  these  equations  and  the  method  of  solution  were  given  in 
Section  4.  2. 

The  numerical  routines  for  solving  Equations  ( 49)  -(52)  have  been  incor- 
porated into  program  SHEATH,  and  a fully  self-consistent  solution  was 
obtained  for  the  parameters  given  in  Figure  13,  The  electron  and  ion 

number  densities  and  drift  velocities  are  shown  as  functions  of  y.  (The 

2 

assumption  that  n+  — const  is  seen  to  be  justified  and  mv^  » kT  over  most 
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of  the  discharge.  ) The  fact  that  ne«  n+  everywhere  indicates  that  the 
anode  is  located  within  the  cathode  fall.  The  situation  is  indicative  of  an 
obstructed  discharge.  If  the  anode  is  moved  back,  the  electron  density 
can  continue  to  rise  until  it  equals  the  positive  ion  density,  at  which  point 
the  negative  glow  begins  to  form. 


As  mentioned  earlier,  the  region  between  the  edge  of  the  cathode  fall  and 
the  head  of  the  positive  column  could  not  be  handled  by  the  code  in  its  original 
form,  because  of  possible  field  reversals.  The  program  has  since  been 
modified  to  allow  a change  of  variables  and  the  use  of  Equation  (55).  This  has 
greatly  increased  the  utility  of  the  code  and  some  typical  results  will  be 
described  in  the  following  section. 


The  Cathode  Fall  and  Negative  Glow 


6.  2 

In  modeling  the  negative  glow,  the  code  was  extended  to  incorporate  a zero  field, 
drift  region.  Electron  loss  mechanisms,  such  as  attachment  and  recombination, 
were  added.  They  do  not  significantly  affect  electron  kinetics  in  the  high  field 
region  of  the  cathode  fall  but  are  crucial  in  the  negative  glow  and  the  transi- 
tion into  the  positive  column.  Likewise,  low  energy  cross  sections,  such 
as  vibrational  and  rotational  excitation  in  N^,  have  little  impact  on  electrons 
within  the  cathode  fall  itself  and  only  come  into  play  at  low  mean  energies. 

One  difficulty  which  came  up  in  the  modeling  concerned  the  presence  of 
field  reversal  at  the  edge  of  the  cathode  fall.  The  electric  field  falls  sharply 
away  from  the  cathode  surface  due  to  a layer  of  nearly  constant  positive  space 
charge.  The  electron  current  density  rises  exponentially  in  this  region  until 
the  electrons  carry  nearly  all  the  current  in  the  negative  glow.  At  the  edge 
of  the  cathode  fall,  the  electron  density  is  rising  much  faster  than  the  current 
density,  since  the  drift  velocity  is  dropping  with  the  field.  Physically  speaking, 
the  secondary  electrons  produced  by  the  propagating  primary  beam  cannot  be 
convected  out  of  the  low  field  region  fast  enough  because  of  their  reduced 
drift  velocity  so  they  pile  up. 

The  difference  between  positive  ion  and  electron  densities  determines  the  slope 
of  the  electric  field  through  Poisson's  equation.  In  the  cathode  fall  region, 
the  positive  ion  density  is  nearly  constant  and  the  electron  density  is  negligible, 
so  the  field  varies  linearly  with  distance.  At  the  edge  of  the  cathode  fall  the 
electron  density  rises  sharply  to  values  which  are  comparable  with  the  positive 
ion  density.  This  causes  the  slope  of  the  electric  field  to  soften  and  perhaps 
to  go  through  zero  and  turn  positive.  The  rate  at  which  electrons  pile  up  at  the 
boundary  between  the  cathode  fall  and  negative  glow  then  determines  the  way 
in  which  the  electric  field  stops  falling  and  joins  with  the  low  field  present  in 
the  negative  glow.  In  particular  the  field  might  easily  go  through  zero  and 
remain  negative  for  a while  before  the  electron  density  comes  up  to  turn  it 
around. 


Since  current  continuity  must  always  be  satisfied,  any  region  of  field  reversal 
must  have  only  a limited  extent.  The  low  energy  electrons  and  positive  ions 
will  follow  the  field  closely,  except  in  cases  of  very  rapid  change.  The  high 
energy  electrons,  then,  are  the  only  thing  available  to  carry  the  current  through 
a field  reversal.  However,  the  current  density  of  high  energy  electrons  is 
falling  monotonically  away  from  the  cathode.  The  balance  is  made  up  by  an 
increase  in  secondary  electron  current.  If  the  field  reverses,  the  secondary 
electron  drift  velocity  will  follow  it  very  closely  resulting  in  a sharp  rise  in 
electron  density  until  it  exceeds  the  ion  density.  Then  the  slope  of  the  electric 
field  goes  through  zero  and  the  field  begins  to  increase.  This  reverses  the 
whole  process  and  the  electron  density  turns  over  and  approaches  the  ion 
density  from  above,  eventually  establishing  a uniform  field  in  the  positive 
column. 

The  existence  of  field  reversal  in  the  negative  glow  presents  serious  computa- 
tional difficulties  in  modelling  the  cathode  region.  In  the  first  place,  the 
localized  potential  well  at  the  edge  of  the  cathode  fall  results  in  trapped  orbitals 
for  low  energy  electrons.  Thus,  the  characteristics  no  longer  terminate  on 
the  electrodes  where  boundary  conditions  are  established  but  form  closed  loops. 
These  orbitals  are  populated  continuously  by  ionization  which  must  be  balanced 
by  recombination  or  some  other  loss  mechamism,  since  no  electrons  can 
escape  by  convection.  But  the  recombination  rate  depends  on  the  positive  ion 
density,  which  is  not  known  until  the  next  step  in  the  iteration  procedure.  In 
addition,  the  presence  of  a large  number  of  low  energy  electrons  requires  a 
much  finer  partitioning  of  the  energy  coordinate  in  this  region. 

In  order  to  obtain  additional  understanding  of  the  electron  kinetics  in  the  cathode 
region  before  addressing  these  computational  questions,  it  was  decided  to  run 
some  typical  cases  using  a field  distribution  which  did  not  go  through  zero. 

This  distribution  consisted  of  two  parts:  a linearly  falling  electric  field  repre- 
senting the  cathode  fall  and  a zero  field  region,  or  drift  region,  representing 
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the  negative  glow.  In  the  drift  region,  we  changed  from  an  energy  to  a position 
coordinate  since  tne  characteristics  are  now  constant  energy  surfaces.  The 
stepsize  is  chosen  as  one -fourth  of  the  minimum  electron  mean  free  path.  The 
cross  sections  are  as  defined  previously,  except  that  we  introduce  aniso- 
tropic scattering.  The  degree  of  elastic  backs cattering  was  chosen  to  decrease 
exponentially  with  increasing  energy. 

The  cross  section  for  dissociative  recombination  of  electrons  with  Ar^,  was 

derived  from  experimental  measurements  of  the  recombination  rate  as  a 
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function  of  electron  temperature.  The  following  energy  dependence  was 
inferred, 

o = 8x10  ^ eV-cm^/f 

The  ionization  and  recombination  rate  constants  in  argon  are  shown  in  Figure  14 
for  an  electric  field  which  falls  from  900  V/cm/Torr  at  the  cathode  to  zero 
at  0.  33  cm  Torr  and  remains  zero  out  to  the  anode  at  0.  83  cm  Torr.  The 
ionization  rate  constant  shows  the  expected  behavior,  starting  from  zero  when 
the  electrons  reach  the  ionization  potential,  then  rising  to  a maximum  and  falling 
with  the  electric  field. 

The  recombination  rate  constant  drops  initially  as  the  primary  group  acceler- 
ates away  from  the  cathode.  As  they  undergo  repeated  inelastic  collisions, 
the  mean  energy  begins  to  decrease  and  consequently  k^  increases.  Just  be- 
fore the  field  goes  to  zero,  ionization  virtually  disappears  and  the  large  group 
of  electrons  near  zero  energy  is  accelerated  by  the  remaining  electric  field. 

As  the  mean  energy  of  this  group  increases,  the  recombination  rate  drops.  In 
the  drift  region,  the  electrons  can  only  lose  energy  and  recombination  once 
again  begins  to  rise.  This  behavior  is  clearly  demonstrated  in  the  evolution 
of  the  energy  distribution  function  to  be  presented  in  the  following. 
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The  Townsend  ionization  coefficient  and  mean  electron  energy  are  shown  in 

Figure  15  for  the  same  conditions  as  above.  The  parameter  a here  is  actually 

the  net  ionization  coefficient  with  recombination  subtracted  out  (fractional 

-7 

ionization  assumed  to  be  n^/N  = 10  ),  Note  that  recombination  becomes 

dominant  in  the  negative  glow  or  drift  region.  The  dashed  curve  represents 
the  equilibrium  ionization  coefficient,  which  is  found  by  letting  the  electrons 
reach  equilibrium  at  the  local  electric  field  strength.  It  has  a maximum  at  the 
cathode  surface  and  goes  to  zero  when  the  field  goes  to  zero.  The  actual 
ionization  coefficient  has  a tail  which  extends  well  into  the  negative  glow,  a 
direct  indication  of  the  nonequilibrium  character  of  the  electron  distribution 
function.  The  mean  energy,  ( , has  the  same  features  as  the  recombination 
rate  but  in  the  inverse  sense. 


The  evolution  of  the  electron  energy  distribution  in  the  cathode  fall  and  drift 
region  is  shown  in  Figure  16.  Note  the  att-STiuation  of  the  primary  group  as  it 
gains  energy  in  the  electric  field.  Just  below  the  primary  group  is  a secondary 
group  which  consists  of  electrons  having  undergone  one  inelastic  collision.  Then 
comes  a long  tail  of  intermediate  energy  secondaries  and  finally  a low  energy 
group.  The  low  energy  secondaries  are  constantly  being  produced  by  the  inter- 
mediate group  until  about  cp  = 140  eV  when  this  group  virtually  disappears.  At 
this  point  the  mean  energy  is  relatively  low  and  there  are  very  few  electrons 
above  ionization  threshold,  but  there  remains  about  a ten  volt  potential  drop 
to  the  edge  of  the  cathode  fall.  Thus,  the  mean  energy  once  again  begins  to 
rise  until  the  field  goes  away  altogether.  This  is  seen  in  Figure  16  as  the  low 
energy  peak  moves  out  to  the  inelastic  threshold.  Those  electrons  which 
cross  the  threshold  and  suffer  inelastic  collisions  begin  to  form  a second 
peak  behind  the  principal  one.  This  process  continues  throughout  the  drift 
region. 

The  electron  current  and  number  density  are  shown  in  Figure  17.  As  before, 
the  current  density  rises  exponentially  in  the  cathode  fall  and  levels  off  in 
the  drift  region.  The  current  actually  falls  slightly  at  the  end  of  the  drift 
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region  due  to  recombination.  The  electron  density  rises  more  rapidly  than 
the  current  density  because  the  drift  velocity  is  falling  with  the  field.  In  the 
zero  field  region  the  drift  velocity  begins  rising  again  because  of  the  presence 
of  an  absorbing  anode.  If  the  anode  were  removed  the  drift  velocity  would  go 
to  zero  exponentially  and,  but  for  the  presence  of  recombination,  the  electron 
density  would  go  to  infinity.  In  reality  this  would  cause  the  electric  field  to 
increase  which  would  reaccelerate  the  electrons,  bringing  their  number  density 
back  down  to  the  ion  density  and  bringing  the  field  to  a constant  value  in  the 
positive  column. 

A similar  case  was  run  in  helium  with  the  electric  field  falling  from  225 

V/cm/Torr  at  the  cathode  to  zero  at  1.  33  cm  Torr  and  the  anode  at  2.  83  cm  Torr. 

This  corresponds  to  the  same  cathode  fall  potential  of  150  V.  The  evolution 

of  the  electron  energy  distribution  is  shown  in  Figure  18.  The  energy  distribution 
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leaving  the  cathode  was  taken  from  the  measurements  of  Hagstrom  and  is 
somewhat  broader  than  that  for  argon.  The  qualitative  features  of  the  distribu- 
tion are  the  same.  The  curve  in  the  lower  right  hand  corner  of  Figure  18  shows 
the  experimental  distribution  of  Gill  and  Webb'’*  taken  near  the  edge  of  the  cathode 
fall.  The  potential  in  their  case  was  about  260  V but  the  general  shape  of  the 
distribution  is  very  similar  to  that  calculated  here. 

With  a complete  knowledge  of  the  energy  distribution  as  a function  of  position 
many  observable  features  of  the  cathode  region  can  be  derived;  for  example, 
the  visible  sidelight  emission.  It  is  assumed  that  this  emission  comes  from 
electronic  levels  of  the  atom  which  are  excited  by  resonant  transitions  from 
the  ground  state.  Using  a general  resonant  type  cross  section,  the  excitation 
rate,  and  hence  visible  intensity  in  an  optically  thin  medium,  was  calculated 
as  a function  of  position.  This  is  shown  in  Figure  19  for  helium.  All  the 
visible  features  in  the  cathode  region  of  a normal  glow  are  evident,  including 
the  thin  cathode  dark  space,  cathode  glow,  and  Crooks  dark  space.  Note  the 
sharp  leading  edge  and  diffuse  trailing  edge  of  the  negative  glow  in  agreement 
with  experimental  observation. 
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Figure  19  The  Cathode  Region  in  Helium 


The  usual  pressure  scaling  in  the  cathode  fall  is  preserved  by  the  nonequilibrium 

formulation  presented  here.  That  is,  the  electric  field  strength  is  proportional 

2 

to  pressure,  p,  the  current  density  scales  as  p and  distances  vary  inversely 
with  p.  Thus,  all  the  characteristics  of  the  cathode  fall  presented  in  Figures  14  to 
19  are  as  valid  at  1 atm  as  at  1 Torr.  The  inclusion  of  attachment  as  an 
electron  loss  mechanism  does  not  destroy  this  scaling.  However,  recombination, 
which  introduces  a term  quadratic  in  the  charged  particle  densities,  does  alter 
the  similarity  laws  in  the  negative  glow.  At  higher  pressures,  where  recom- 
bination is  more  significant,  the  length  of  the  negative  glow  is  expected  to  be 
reduced  from  the  usual  scaling  prediction. 


7.0 


SHEATH  CHARACTERISTICS 


An  investigation  of  the  voltage-current  characteristics  of  the  cathode  sheath 
using  the  computer  program  SHEATH  has  shown  that  there  is  a range  of 
parameters  where  the  gross  electrical  features  are  approximated  by  the 
Townsend  model  of  ionization  with  a linearly  varying  electric  field.  Devia- 
tions occur  at  very  high  and  very  low  values  of  the  cathode  fall  potential, 

V . As  V increases  for  a given  cathode  fall  distance,  d , the  electrons 
c c c 

are  accelerated  quickly  to  energies  above  the  peak  of  the  ionization  cross 
section  and  current  multiplication  in  the  sheath  is  reduced.  This  behavior 
is  seen  in  the  numerical  solution  but  not  in  the  equilibrium  model.  At  low 
values  of  , on  the  order  of  several  times  the  ionization  potential,  the 
free-fall  region  near  the  cathode,  where  no  ionization  takes  place,  becomes 
important.  This  region  is  not  considered  in  most  equilibrium  models. 

At  intermediate  values  of  V , near  the  normal  operating  conditions, 
a simple  equilibrium  model  with  linearly  varying  field  has  been  used  to 
estimate  gas  heating  in  the  cathode  sheath.  Temperature  profiles  are  calcu- 
lated and  it  is  found  that  thermal  conduction  to  the  cathode  is  more  important 
than  convective  cooling  for  moderate  Mach  number  flows. 

7.  1 Linear  Field  Approximation 

It  has  been  observed  experimentally  and  confirmed  in  numerical  simulations 
of  the  cathode  fall  region  that  the  electric  field,  and  consequently  the  dis- 
charge power  density,  vary  linearly  with  distance  from  the  cathode.  This 
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is  true  for  both  equilibrium  and  nonequilibrium  models  with  constant  gas 
density  as  well  as  for  a model  incorporating  the  full  gas  dynamic  equations 
of  a compressible  boundary  layer.  * Using  this  empirical  observation  as  a 
starting  point,  it  is  possible  to  construct  a very  simple  yet  accurate  picture 
of  the  cathode  fall  including  gas  flow. 


First  we  will  derive  closed-form  expressions  for  the  cathode  fall  distance 
and  potential  drop  as  a function  of  the  current  and  gas  pressure.  Then,  the 
power  density  will  be  calculated  to  determine  the  temperature  rise  in  the 
thermal  boundary  layer  under  steady-state  conditions.  This  temperature 
rise  corresponds  to  a gas  density  defect  at  the  entrance  to  the  positive 
column  which  can  lead  to  thermal  instabilities  in  the  plasma.4* 

The  electron  current  density  in  the  cathode  fall  satisfies  the  following  con- 
tinuity equation, 


— a je 


(71) 


where  a is  Townsend's  first  ionization  coefficient  and  y is  the  distance  from 
the  cathode.  The  processes  of  electron  loss,  such  as  attachment  and  re- 
combination, are  negligible  compared  to  ionization  in  the  cathode  fall.  The 
solution  of  Equation  (71)  is 

je  = jeoeXPk7a(y)dy]  (72) 

where  j is  the  electron  current  density  leaving  the  cathode. 


In  Townsend's  original  experiments,  electrons  produced  by  illumination  of 
the  cathode  were  allowed  to  drift  in  a uniform  electric  field  and  the  current 
growth  was  measured  as  a function  of  distance.  If  the  electrode  separation 
is  not  too  small,  the  electrons  spend  most  of  their  time  in  equilibrium  between 
the  field,  E,  and  the  gas  density,  N,  and  a depends  only  on  the  ratio  E/N. 

This  dependence  has  the  form 

a = Apo  exp(-  BpQ/E)  (73) 

where  PQ  is  the  pressure  reduced  to  0°C  and  is  proportional  to  N.  A and  B 
are  constants. 
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In  the  cathode  fall  region,  the  electric  field  and  gas  density  are  changing 
so  rapidly  with  position  that  the  electrons  never  reach  equilibrium,  and 
it  is  not  possible  to  relate  a to  any  local  variable.  (It  has  been  suggested^ 
that  the  mean  energy,  e , and  through  it  a , may  be  expressed  as  functions 
of  the  average  field,  V(y)/y.  This  approach  was  later  modified^to  account 
for  the  free-flight  region  near  the  cathode.  ) However,  the  coefficient  a 
appears  in  Equation  (72)  under  an  integral  over  position.  In  the  present  study, 
it  has  been  found  that  the  use  of  Equation  (73)  in  (72)  with  a linear  field  de- 
pendence gives  an  accurate  prediction  of  the  electron  current  growth  in  the 
cathode  fall.  In  other  words,  the  expression  for  a (E/p  ) may  not  apply 
locally,  but  in  an  integral  sense,  over  the  whole  region,  it  has  some  validity. 


If  the  electric  field  is  given  by  the  expression, 

E,y)  = Eo(*  - ar)  <?4> 

where  Eq  is  the  electric  field  at  the  cathode  and  d is  defined  as  the  cathode 

c 

fall  distance,  then 

-S/E°- !-/,(#)]  <75> 

where  Ej  is  the  exponential  integral  function,  A = Ap  and  B = B p . 

° o* 

The  electron  and  positive  ion  current  densities  at  the  cathode  are  related 
by 

jeo  = yj+  o (76) 

where  yis  the  secondary  electron  emission  coefficient.  At  the  edge  of  the 
cathode  fall,  the  positive  ion  current  density  is  negligible  so 
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j a(y)dy  = Ad 


and  In  fi  (d  )/i  ] = In  (1  + 1 17).  The  cathode  fall  distance  can  therefore 

e c eo 


be  written  from  Equation  (75)  as 


ln(  1 + 1 /Y) 
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(78) 


The  current  density  can  be  found  from  Poisson's  equation. 
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evaluated  at  the  cathode  surface 
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where  v+  is  the  positive  ion  drift  velocity. 


In  the  cathode  fall  region,  the  positive  ion  mobility  is  limited  by  charge 
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exchange.  If  the  cross  section  for  this  process  is  nearly  constant  then 


v+  = k 


• (t) 


(81) 


If  we  now  define 


&). 


(82) 


then  from  Equation  (80) 
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(83) 


This  expression  is  valid  for  any  field  dependence,  not  just  the  linear  depen- 


dence of  Equation  (74).  Equation  (83)  agrees  very  well  with  experimental 


observations  in  helium  at  1 Torr  and  with  numerical  predictions.' 
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We  now  have  the  current  density  and  cathode  fall  distance  expressed  in  terms 

of  the  physical  constants,  A,  B,  and  7,  and  the  parameter,  E . With  the 

o 

cathode  fall  potential  given  by  V = 1/2  E d , we  can  plot  p d and  V as 
2 c oc  o c c 

functions  of  j/pQ  (Figure  20).  The  minimum  in  the  cathode  fall  potential  is 

defined  as  the  normal  operating  condition,  with  the  properties 


Jn 

0.  57  AB3/2  € k 
o + 

1 + Y 

(84) 

Po2 

ln(l  + 1/7) 

V 

n 

3 ln(l  + l/Y) 

(85) 

P d = 
o n 

“■  In  (1  + 1/7) 

(86) 

These  results  for  constant  temperature  are  presented  in  many  standard  texts 
on  gaseous  electronics. 

For  a discharge  in  argon  with  an  iron  cathode,  7 = 0.04,  A = 14  cm’1  Torr’1, 
B = 180  V cm"1  Torr"1,  and  k + = 8.25  x 103  cm3/2  Torr1/2  V_1/2  s'1. 

In  this  case 

~ - 4.  5 x 10  A cm  Torr 

Po 

V = 126  V 

n 

p d = 0.  86  cm  Torr 

o n 

These  values  are  in  close  agreement  with  the  numerical  calculations  of  Ward3® 
even  though  he  used  a different  expression  for  a (E/p  ).  However  the 
normal  current  density  is  considerably  low  and  the  distance  somewhat  high 
in  comparison  with  experiments.  (The  experimental  values  can  be  obtained 
with  a different  choice  of  physical  constants.  ) 


- — .•  M • . 
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7.  2 Gas  Heating  in  the  Cathode  Sheath 

The  reduced  pressure  appearing  in  Equation  (83)  is  measured  at  the 
cathode  surface.  Because  of  gas  heating  in  the  cathode  fall  region,  Pq  may 
be  considerably  different  from  the  pressure  in  the  positive  column.  The 
temperature  rise  in  the  cathode  fall  can  be  determined  from  a knowledge  of 
the  power  distribution  and  the  thermal  conductivity  of  the  gas.  Since  the 
current  density  is  independent  of  position,  the  discharge  power  density,  is 
given  by 


-Po{ 


'-f- 

C 


where  P = jE  . Essentially  all  of  the  power  dissipated  in  the  cathode  fall 
goes  into  heating  the  positive  ions  which  in  turn  heat  the  gas  through  charge 
transfer  collisions.  Thus,  Equation  (87)  can  be  considered  as  the  rate  of 
heat  addition  to  the  gas. 

In  a discharge  channel  with  flow  parallel  to  the  cathode  surface,  a boundary 
layer  is  established  due  to  the  finite  viscosity  of  the  medium.  The  gas 
velocity  in  the  boundary  layer  goes  from  zero  at  the  cathode  to  a maximum 
value  in  the  free  stream.  For  gas  flows  near  standard  density  • Pfi,  the 
thickness  of  this  layer,  6 , is  much  greater  than  the  normal  cathode  fall 
distance,  and  increases  in  the  downstream  direction.  For  example,  in  the 
case  described  in  Reference  ( 1)  with  a free  stream  density,  pg  = 0.  3pg, 
the  boundary  layer  displacement  thickness  was  0.2  cm  at  a distance  10  cm 
downstream  of  the  nozzle.  The  cathode  fall  thickness  was  found  to  be  0.007  cm. 
Thus,  the  cathode  fall  is  well  within  the  laminar  sublayer,  and  for  subsonic 
flows  at  least,  the  effects  of  convective  heat  transfer  in  this  region  may  be 
neglected. 

We  now  define  the  heat  conduction  problem  by  considering  a solid  cathode 
of  thickness  t interfacing  the  gas  at  y = 0.  The  temperature  of  the  back 
face  of  the  cathode  has  the  constant  value  at  y = -t.  Heat  addition  to 
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the  gas  is  defined  by  Equation  (87)  and  the  heat  flux  into  the  stream  at 
y = d^  is  represented  by  q^.  The  thermal  conductivity  of  the  gas  is  assumed 
to  have  a power  law  dependence  on  temperature,  i.  e. 


K 

g 


(88) 


where  is  some  reference  temperature,  here  taken  to  be  the  temperature 
at  y = -t.  The  gas  temperature  in  the  cathode  fall  satisfies  the  equation. 


a / st\  • 

Tf  + QW  = 0 <89) 

where  Q(y)  = P_(y).  In  the  cathode  itself 
a 

d2  T 

K = 0 (90) 

dy 

where  k c , the  thermal  conductivity  of  the  cathode  material,  is  assumed 
constant.  The  temperature  distribution  is  then, 


K (T  - T ) = ( — Q d - q W + t) 
c o V 2 ° c c' 


(91) 


for  -t  < y < 0 and 


K T 
o o 

1 + a 


(92) 


for  0 < y <dc»  where 

*c|Tw-Tol=  (4-  Vc-O1  • «”» 
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The  temperature  at  the  edge  of  the  cathode  fall,  T , is  specified  by  Equa- 


tion (92)  with  y = d . The  heat  flux  into  the  free  stream  is  given  by  the 


usual  expressions  for  local  heat  flux  in  laminar  or  turbulent  flows  over 
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heated  flat  plates,  i.  e. 


Laminar:  q (x)  = 0.  332  k P ^ R ^ (T  - T ) — — 

c x c a x 


(94a) 


Turbulent:  q (x)  = 0.0288  k R 4/5  (T  - T ) — 

c x c a x 


(94b) 


Here  the  thermal  conductivity,  k,  was  assumed  constant,  P is  the  dimension- 
less Prandtl  number,  R the  Reynold's  number  (Ux/v),  and  T = T + \/p  ^ — 

x a ® 2 c 

is  the  adiabatic  wall  temperature.  In  both  cases,  the  heat  flux  increases  ^ 


with  T and  decreases  with  x,  the  distance  measured  from  the  leading  edge. 

C 5 -1/2 
However,  for  laminar  flow  (Rx<  5 x 10  ) q^(x)  , while  for  turbulent 

flow  (R  > 5 x 10^)  q (x)  — x . 
x c 


An  estimate  will  now  be  made  of  the  relative  importance  of  q in  determining 


the  temperature  distribution  in  the  cathode  fall.  For  a turbulent  flow  with 


a = 1 and  T = T , the  temperature  at  y = d is  found  from  Equation  (22)  to 
a o c 


ac  (T  - T HI  + 0.  0058  R 

O C O \ X 


4/5  dc 


4“  Q d 2 

6 o c 


The  temperature  without  flow  is  obtained  by  setting  Rx  = 0.  Thus  convective 

cooling  can  be  neglected  when  U « U = — / X \ ^ ^ . When  U = U.,.  , 

dc  \ dc  / 

the  temperature  at  the  edge  of  the  cathode  fall  is  reduced  by  a factor  of  two. 


The  critical  Mach  number,  M.,,  , is  given  by 


c p d 
“ w c 


(t) 


where  p and  c^  are  the  free  stream  viscosity  and  sound  speed,  respectively, 
anc*  = P (P  /760  Torr)  is  the  gas  density  at  the  wall. 

w SO 


77 


For  nitrogen  at  T = 300  K this  condition  becomes 

5/4 


M 


/ y_  i 2 \ ’ 

+ ^ ) = 34.5  (£*) 


1/4 


(97) 


with  p in  amagat  and  x in  cm.  The  critical  Mach  number  at  x = 10  cm 

00 

in  at  standard  density  is  M.,.  = 5.  6. 

The  temperature  distribution  defined  by  Equation  (92)  agrees  very  well  with 
the  numerical  results  of  Reference  (1).  When  convective  cooling  is  ignored 
the  temperature  rise  in  the  cathode  fall  region  is  from  Equation  (92). 


K T 

0 o 

1 + a 


X 

T 


J +a 


T \ 1+a 

w I 


-j—  Q d 2 

6 o c 


(98) 


Here  we  see  the  importance  of  cooling  the  cathode  in  order  to  minimize  the 
density  defect  caused  by  gas  heating.  However,  cooling  the  cathode  to  a 
temperature  much  below  T can  have  only  a limited  effect  if  AT  » T . 


The  boundary  layer  on  a flat  plate  makes  a transition  from  laminar  to  turbu- 

46 

lent  flow  at  a Reynold's  number, 


U x 


-3  x 10' 


(99) 


where  U is  the  free -stream  velocity,  x is  the  distance  from  the  leading 
00 

edge  of  the  plate  and  v is  the  kinematic  viscosity.  In  nitrogen  at  STP 
2 

V - 0.  133  cm  /sec,  and  the  transition  occurs  at  about 

1 . 2 cm 


crit 


M 


(100) 


shere  M U /c  is  the  free -stream  Mach  number.  For  x > x ..  momen- 
ts ao  crit, 

• and  heat  transfer  in  the  boundary  layer  are  dominated  by  turbulent 

, 46 

. ..  I he  thickness  of  the  layer  is 


= 0.  37  x R 


-0.2 


(101) 
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corresponding  to  the  point  at  which  the  longitudinal  velocity  reaches  99  per- 

0 8 

cent  of  its  free -stream  value.  For  a given  velocity,  5 grows  as  x * until 
the  entire  flow  channel  is  turbulent. 


The  normal  thickness  of  the  cathode  dark  space  in  a glow  discharge  is 

• ,,  52 

typically 


d p 
n o 


0.3  - 1.3  cm  Tor r 


(102) 


where  p is  the  reduced  pressure  at  the  cathode  surface,  i.e., 

273  K 


(103) 


In  a nitrogen  discharge  at  STP  with  an  iron  cathode 

0.  42  cm  Torr 


52 


d p 
n o 


(104) 


and 


-4 

5.  5 x 10  cm. 


On  the  other  hand,  we  see  from  Equation  (101)  that  « > 0.036  cm/M  in  the 
turbulent  regime.  Thus,  under  normal  conditions  at  atmospheric  pressure, 
the  cathode  sheath  is  embedded  well  within  the  turbulent  boundary  layer. 


For  a given  temperature  the  kinematic  viscosity  is  inversely  proportional  to 

pressure,  so  5 varies  as  p On  the  other  hand,  d varies  as  p \ Thus, 

n 

at  low  pressure,  the  ratio  of  d to  * is  small  only  at  very  large  x,  but  in 
this  case  the  transition  to  turbulence  also  occurs  at  larger  values  of  x. 


Written  ii.  terms  of  X . we  have 

crit 


d 

n 


(105) 


where  the  constant  is  independent  of  pressure. 


In  the  region  very  close  to  the  wall,  the  velocities  are  so  small  that  viscous 

forces  continue  to  dominate  over  inertia  forces.  This  is  known  as  the  laminar 
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sublayer.  The  thickness  of  this  layer  on  a flat  plate  is 


= 29.  5 x R 


(106) 


For  a given  velocity,  5 grows  asx  ' , or  much  more  slowly  than  a.  Using 

^ > -4 

Equation  (100)  in  Equation  (106)  we  see  that  6^  ~ 4.  2 x 10  cm/M  in  nitrogen. 

Thus,  for  subsonic  flows,  the  normal  cathode  dark  space  and  laminar  sublayer 

-0.  9 

have  roughly  the  same  extent.  The  pressure  scaling  of  * is  as  p , simi- 
lar to  d . Written  in  terms  of  x for  nitrogen  at  0°C,  we  have 
n crit 


— = 1 


•■ft*) 


(107) 


where  once  again  the  constant  is  independent  of  pressure.  For  argon  at  0"C, 
it  becomes  1. 2,  while  for  helium  at  the  same  temperature,  it  is  0.  5. 

In  the  laminar  region,  the  transport  of  heat  normal  to  the  flow  takes  place  by 
conduction.  This  can  be  seen  by  introducing  the  dimensionless  form  of  the 
energy  equation  with  constant  fluid  properties 


»h°  :'h0  1 l^r  , , E • , n 

“ — +V-T  = PR  [ 2 + 2 ] + R {+Q 

J 7 x \ >x  }y  / x 


(108) 


2 Vc 

where  $ is  the  viscous  dissipation  function,  P - n is  ll>e  Prandtl 

5Y  z 

number,  and  E = — , is  the  Eckert  number. 

c (iT) 

p o 

In  the  definitions  above,  c is  the  heat  capacity  per  unit  mass,  K is  the  thermal 

P 

conductivity,  and  („T)  is  the  difference  between  the  wall  temperature  and  the 
maximum  temperature  in  the  boundary  layer.  The  lengths  in  Equation  (108) 
are  normalized  to  x and  the  velocities  to  U^.  The  rate  of  heat  addition  is 
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normalized  as  follows- 


Q x 


Pc  U (AT) 
p “ o 


(109) 


The  orders  of  magnitude  of  the  normalized  velocities  in  Equation  (108)  are 
found  from  the  solution  of  the  laminar  boundary  layer  equations  to  be 


u 


and 


The  derivatives  of  h and  T with  respect  to  x are  of  order  1.  The  derivatives 
with  respect  to  y are  estimated  as  follows:  Since  the  heat  addition  takes  place 
within  a sheath  of  thickness,  d,  the  temperature  (and  h^)  will  reach  a maxi- 


mum value  at  y , such  that  0 < y < d.  Thus,  in  the  region  0 <■  y < y 


max 


and 


ay 

^2 
2 


max 

x 
d 


max 


5y 

Later  we  will  see  that  y^  ~d  and  these  approximations  are  valid  through- 
out the  sheath. 

Now,  replacing  the  terms  in  Equation  (108)  with  their  orders  of  magnitude, 
we  have 


1 + 


& _ 1 / x \ _E_  x" 

d ' PR  ( 1 + 2 ) + R 2 

x \ d / x 


+ Q 


(110) 


The  terms  on  the  left  in  Equation  (110)  are  both  of  order  unity  for  the  normal 
cathode  sheath.  The  first  term  on  the  right  is  small  compared  to  the  second, 
i.e.,  there  is  little  or  no  heat  flux  along  the  surface.  Using  Equations  (101) 
and  (105),  the  second  term  becomes 


“t  (rH 

PM  \ crit  f 


PM 

where  P ~ 1 for  gases.  So  for  subsonic  flow  and  x >x  the  convective 

crit 
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terms  on  the  left  of  Equation  (108)  are  negligible  compared  to  the  conductive 
term  on  the  right.  The  ratio  of  viscous  heating  to  thermal  conduction  is  of 
the  order  E,  which  is  small  when 


U 


T) 


0 


(111) 


Equating  the  remaining  two  terms  in  Equation  (110)  we  find,  for  P = 1, 


(uTP 


Qd 


(112) 


In  the  cathode  sheath,  the  principal  heat  addition  is  by  collisions  with  posi- 
tive ions,  which  account  for  most  of  the  current  in  a cold  cathode  sheath. 
Thus, 

Q ~ jE  (113) 

where  j is  the  current  density  and  E the  electric  field.  Still  dealing  in  orders 
of  magnitude,  we  can  rewrite  Equation  (112)  as 

WT»o~i¥  <114> 


where  V r Ed  is  the  drop  in  potential  across  the  sheath.  Typical  values  of 
(AT)  for  various  gases  on  an  iron  cathode  are  shown  in  Table  I at  1 Torr 
and  30CPK. 


Table  I 


He 

H2 

Ar 

N2 

(AT)0 

0.  26 

8.  7 

18.  5 

139. 

(CK) 

These  were  obtained  by  using  normal  values  of  current  density,  potential  and 
sheath  thickness. 

Now  that  we  have  an  equation  describing  heat  dissipation  in  the  sheath,  we 
can  construct  a model  which  couples  the  field,  the  plasma,  and  the  gas  in  a 
self-consistent  manner.  In  order  to  examine  the  effects  of  gas  density 


variations  within  the  sheath  on  the  sheath  characteristics,  we  begin  with  a 
simple  model  which  treats  the  electrons  as  being  in  equilibrium  with  the 

' 

field.  It  has  been  shown  in  Section  6 that  this  is  not  a good  approximation 
in  treating  the  fine  details  of  sheath  structure  such  as  electron  density  oscil- 
lations near  the  cathode.  However,  for  determining  the  electric  field  vari- 
ation with  position  and  the  rate  of  heat  addition  to  the  gas,  it  is  quite  an 
adequate  approximation  in  the  vicinity  of  the  normal  operating  point. 

The  equations  to  be  solved  for  the  steady  state  are  simply 

t — — = S+aj-rnn  (115) 

}y  Je  e e + 


j = e n v 

J c e e 
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The  various  constants  for  nitrogen  are 


A 
B 

H P 

e o 

k 

+ 

(T) 
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e 

The  numerical  integration  is  started  at  the  edge  of  the  sheath,  where  n+ 

and  proceeds  until  the  boundary  conditions  j = and  T = T^  are  satisfied 
at  the  cathode  surface.  Since  the  two  conditions  must  be  satisfied  simulta- 
neously, a "shooting " technique  is  employed,  where  the  initial  temperature 
is  varied  until  convergence  is  achieved. 

With  regard  to  pressure  scaling,  we  make  the  following  observations.  In 
the  absence  of  recombination  and  gas  heating,  the  similarity  variables  (i.e., 
variables  which  make  the  system  of  equations  independent  of  pressure)  are 
S/p  , j/p  , E/p,  and  yp.  Using  these  variables,  the  recombination  term  in 
Equation  (115)  and  the  heat  addition  term  in  Equation  (117)  are  directly  pro- 
portional to  pressure.  Recombination  is  only  significant  in  the  negative  glow 
and  positive  column  and  does  not  affect  the  sheath  characteristics  appreciably. 
However,  as  pointed  out  in  the  last  section,  gas  heating  is  very  important 
even  at  low  pressures.  Thus,  the  usual  pressure  scaling  no  longer  holds. 

If  K were  constant,  then  temperature  changes  in  the  cathode  sheath  would 

scale  directly  with  pressure.  Since  K increases  with  T,  the  actual  scaling 

should  be  less  than  linear  with  pressure.  However,  the  heating  term  is 

modified  bv  the  variation  of  p in  the  sheath  and  it  will  be  shown  below  that 
' o 

the  effects  of  variable  K and  variable  pQ  tend  to  offset  one  another,  such 
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that  AT  still  increases  directly  with  pressure.  This  fact  represents  the 
major  stumbling  block  to  operating  cold  cathode  discharges  at  high  pressure. 

The  electric  field  in  the  sheath  is  plotted  in  Figure  21  as  a function  of  position 

for  three  different  pressures.  Scaling  variables  are  used  so  that,  if  p were 

o 

constant,  all  three  curves  would  coincide.  The  deviations  are  minor  for  pres- 
sures up  to  300  Torr.  Note  in  particular  that  E/p  at  the  cathode  surface  is 
practically  constant.  The  cathode  fall  potential  is  only  slightly  higher  at 
300  Torr  than  at  10  Torr,  but  increases  significantly  at  1 atm. 


The  temperature  distribution  in  the  sheath  is  shown  in  Figure  22.  The  flow 
velocity  is  assumed  to  be  small  so  that 
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By 
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(118) 


This  is  also  true  under  static  conditions  if  the  electrode  spacing  is  large 
compared  to  the  sheath  dimension.  Notice  that  the  width  of  the  thermal 
boundary  layer  is  less  than  the  cathode  fall  distance  and  that  the  tempera- 
ture reaches  its  maximum  value  well  within  the  sheath.  At  300  Torr,  the 
temperature  at  the  plasma -sheath  boundary  is  950  K and  the  gas  density  is 
only  one  third  of  its  free  stream  value.  In  an  e-beam  sustained  device 
where  the  electric  field  is  determined  by  the  applied  voltage,  the  E/p  or 
E/N  at  the  edge  of  the  sheath  will  be  three  times  as  high  as  in  the  core  flow. 
Thus,  if  E/N  must  be  maintained  below  a critical  value  where  volume  ion- 
ization begins  to  compete  with  the  e -beam  source,  the  temperature  in  the 
boundary  layer  presents  a severe  limitation  on  the  applied  voltage,  and  hence 
the  power  density,  which  can  be  achieved  before  instabilities  occur. 


In  the  transient  case,  where  a voltage  pulse  is  applied  at  t = 0,  the  tempera- 
ture distributions  of  Figure  22  will  take  some  time  to  develop.  There  are 
two  time  scales  of  importance  here:  The  time  required  to  heat  the  gas  in 
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a: 


the  sheath,  T , and  the  time  it  takes  a pressure  disturbance  to  propagate 
across  the  sheath,  T2.  The  first  characteristic  time  is  approximately 


jE 


(AT) 


(119) 


where  c is  the  specific  heat  per  unit  mass  at  constant  volume.  Taking  for 
v r 


nitrogen  pc  =3  Nk  and  using  Equation  (114),  we  have 


1 


3_Nk  d2 


(120) 


where  N is  the  initial  gas  density.  In  nitrogen  at  STP  and  normal  current 
density  on  an  iron  cathode, 

t ~ 0.65  p^ec 

This  value  scales  as  p \ so  that  a 1 Torr  ^ 0.  5 msec. 

The  second  characteristic  time  is  just 

d 


(121) 


Again  in  nitrogen  at  STP,  T ~ 0.  017  Msec.  This  also  scales  as  p . Since 


the  pressure  in  the  sheath  equilibrates  on  a time  scale  corresponding  to  t^, 
which  is  much  less  than  T , the  assumption  of  constant  pressure  applies  to 
some  degree  in  the  transient  regime  as  well  as  in  steady  state. 


It  should  be  pointed  out  here  that  if  the  discharge  is  operated  at  current  den- 
sities greater  than  j , then  d is  less  than  d , and  both  and  t will  be  smaller 

n n i <- 

than  their  normal  values.  Attempts  to  operate  at  current  densities  less  than 

j will  result  in  bunching  of  the  current  at  the  cathode  such  that  the  local 
n 

i = i . This  is  a two-dimensional  effect  which  has  not  been  treated  under  the 
n 

present  contract  due  to  time  limitations. 

The  temperature  rise  across  the  sheath  is  plotted  in  Figure  23  as  a function 
of  pressure.  Note  the  nearly  linear  dependence  as  suggested  earlier.  The 
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appear!  ice  of  these  temperature  increments  and  associated  density  defects 
at  the  edge  of  the  cathode  sheath  on  time  scales  as  defined  above  must  have 
serious  consequences  for  the  stability  of  discharge  devices  at  high  pressure. 
Although  the  adiabatic  propagation  of  density  disturbances  at  the  local  sound 
speed  is  too  slow  to  account  for  many  of  the  observed  instabilities  in  TEA 
discharges  and  e-beam  sustained  devices,  these  disturbances  may  be  driven 
by  locally  enhanced  ionization  and  gas  heating  at  much  greater  velocities. 
Therefore,  the  coupled  field,  plasma,  and  gas  equations  must  be  treated  self- 
consistently  and  time  -dependently  in  order  to  deal  with  this  situation.  In 
addition  to  longitudinal  instabilities,  there  are  other  instabilities  resulting 
from  transverse  nonuniformities  in  the  electric  field  or  external  ionization 
which  cannot  be  dealt  with  under  the  present  framework  but  which  must  be 
considered  in  a general  stability  analysis. 

The  full  temperature  distribution  in  the  cathode  sheath  and  boundary  layer  is 
shown  in  Figure  24  for  three  different  Mach  numbers.  The  sheath  dimension 
has  been  expanded  relative  to  the  boundary  layer  thickness  in  order  not  to 
obscure  the  detail  in  this  region.  The  temperature  and  heat  flux  at  the  sheath - 
boundary  layer  interface  are  continuous.  Note  the  minimal  cooling  provided 
by  low  Mach  number  flows.  The  relative  density  change  across  the  boundary 
is  actually  greater  at  higher  Mach  numbers  due  to  the  adiabatic  cooling  of 
the  core  flow. 
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8.  0 STABILITY  ANALYSIS 


T 


The  temperature  and  electric  field  at  the  cathode  surface  are  critical  in 
determining  the  stability  of  the  discharge.  If  thermally  enhanced  electron 
emission  becomes  significant,  the  uniform  glow  collapses  into  an  arc.  This 
occurrence  can  be  forestalled  by  proper  cooling  and  conditioning  of  the 
cathode.  For  a given  pressure  and  electrode  configuration,  however,  there 
exists  a maximum  current  density  for  stable  operation.  Electrode  cooling 
and  thermal  stability  are  covered  in  Section  8.  1. 

The  effects  of  field  tailoring  and  auxiliary  ionization  on  sheath  characteris- 
tics and  discharge  stability  are  examined  in  Section  8.  2.  It  is  found  that 
field  tailoring  in  e -beam  sustained  devices  is  effective  in  countering  thermal 
and  chemical  changes  which  occur  in  the  flow  direction.  These  changes 
occur  principally  in  the  positive  column,  but  they  can  modify  the  current 
density  distribution  which  affects  the  sheath.  Auxiliary  ionization  typically 
has  a negligible  effect  on  sheath  characteristics,  but  it  can  lead  to  thermal 
instabilities  at  the  edge  of  the  cathode  fall  which  propagate  into  the  positive 
column. 

8.  1 Electrode  Cooling  and  Thermal  Stability 

In  the  last  section,  it  was  assumed  that  the  cathode  surface  was  maintained 
at  a temperature  T = 300  K.  In  practice,  the  cathode  must  be  cooled  on 
the  side  away  from  the  discharge,  for  example,  by  passing  a fluid  with  high 
heat  capacity  over  the  surface.  The  temperature  at  the  cathode-gas  inter- 
face then  depends  upon  the  thermal  flux  incident  from  the  discharge.  The 
effect  of  raising  the  temperature  at  the  cathode  surface  is  to  decrease  the 
density  of  adjacent  gas  and,  if  the  current  density  is  held  constant,  to  in- 
crease the  thermal  flux.  This  process  leads  to  an  instability  in  which  the 
cathode  temperature  and  electric  field  at  the  surface  increase  without  bound. 
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The  combination  of  increasing  temperature  and  electric  field  eventually 
changes  the  electron  emission  process  from  one  which  is  independent  of  T 
to  one  which  increases  with  T,  e.  g. , thermionic  field  emission.  When  the 
emission  of  electrons  is  an  increasing  function  of  T,  the  slope  of  the  voltage  - 
current  characteristic  for  the  sheath  goes  negative  and  the  discharge  col- 
lapses into  an  arc.  In  the  following  paragraphs  we  will  show  that  this  process 
is  inevitable  once  the  local  current  density  is  raised  above  a critical  value 
which  depends  on  pressure. 


At  low  temperature,  the  principal  mechanisms  for  emission  of  electrons 
from  the  cathode  are  ion  impact  and  photoemission.  Under  ion  impact,  the 
electron  current  at  the  cathode  is  related  to  the  ion  current  by 


e,  1 


= v+ 


(122) 


where  y.  is  the  electron  yield  per  ion.  The  factor  is  a function  of  ion 
species  and  energy  as  well  as  surface  material  and  condition.  In  the  high 
field  region  of  the  cathode  sheath  the  ion  species  is  typically  the  singly  ion- 
ized component  of  the  parent  gas  with  an  average  energy 

e.  ~ E/NQ  (123) 

i 

where  Q is  the  cross  section  for  symmetric  charge  exchange.  In  argon  with 

E /p  ~ 1000  V cm"'  Torr'',  e.  ~ 5 eV  and  in  helium  with  E /p  ~230  V cm 
o ro  1 o o 

Torr‘1,  ~ 2.  5 eV.  The  yields  for  these  ions  on  atomically  clean  tungsten 
are^  y.  (Ar+)  = 0.  10  and  y^  (He+)  = 0.29.  On  tungsten  covered  with  a mono- 
layer  of  N2>  yi  (Ar+)  = 0.  035  and  (He+)  = 0.  17.  Thus  the  treatment  of  the 
surface  is  critical  in  determining  accurate  sheath  characteristics. 


In  addition  to  ions  incident  on  the  cathode  surface  there  may  also  be  photons 

with  sufficient  energy  to  eject  an  electron.  The  threshold  frequency  for  photo - 

emission,  v i is  defined  by 
o 

= CD  (124) 
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where  cp  is  the  work  function  of  the  material.  This  ranges  from  2 eV  to  5 eV 

for  most  metals.  An  expression  for  the  photoelectric  current  was  derived 
54 

by  Fowler 


e,  ph 


= CT 


f \Hv-v  ) 

] i oox  hh 


(125) 


where  C is  a constant  that  depends  on  the  properties  of  the  solid,  I is  the 
incident  radiation  intensity,  and  X is  a universal  function  which  Fowler 
expressed  in  series  form. 

The  intensity,  I(p),  is  dependent  on  the  gas  and  the  operating  conditions  and 
geometry  of  the  discharge.  In  general,  we  can  write 


I(l/)  = "n(p,  E)  j EX(p) 


(126) 


where  Ti  is  the  efficiency  for  producing  photons  of  frequency  u in  the  discharge 
and  l is  a characteristic  length.  If  the  radiation  is  not  trapped  by  the  gas, 
then  X is  defined  by  the  geometry  of  the  discharge  region.  If  there  is  signifi- 
cant trapping,  then  L is  related  to  the  trapping  distance. 

For  a given  gas  and  applied  electric  field,  the  photoelectric,  current  is 
directly  proportional  to  j,  i.e., 


e,  ph 


= 'Vh  j 


(127) 


where  y is  a function  of  E and  various  geometrical  factors.  Now  we  can 
ph 

relate  the  total  electron  current  to  the  ion  current  at  the  cathode  surface  by 


and  define  an  effective  V 


vi  + Vph  J 


1 - v 


^i + Yph 


(128) 


vh  ich  includes  both  ion  and  photon  processes.  In  a gas  where  the  uv  emission 
comes  principally  from  transitions  between  excited  or  ionic  states,  there  is 
little  trapping  and  V may  be  much  greater  than  y^. 
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The  cathode  fall  potential  at  a given  current  density,  and  hence  the  heat 
loading  in  the  sheath,  is  a sensitive  function  of  V.  In  particular,  the  normal 
current  density  increases  with  V.  Therefore  the  choice  of  gas  and  cathode 
material  and  the  treatment  of  the  cathode  surface  can  have  a significant 
effect  on  the  operating  characteristics  of  the  discharge. 


(129) 


The  heat  flux  to  the  surface  of  the  cathode  is  just 

K ^ = jV 

ay  y=o  J 

2 2 

For  a given  Y the  cathode  fall  potential,  V,  is  a function  of  P = iT  / p , 

w w 

where  T is  the  temperature  at  the  cathode  surface.  If  the  cathode  thick- 
w 

ness  is  t and  the  back  face  is  held  to  a temperature,  T , then  T^  is  given 
by 


— (T  - T ) = jV(§  ) 
t w o w' 


(130) 


where  is  the  thermal  conductivity  of  the  cathode.  Solving  Equation  (30) 
for  j in  terms  of  the  similarity  parameter  §,  we  have 

1/2  _ 1/2 


.3/: 

J 


w 


- 


m) 


(131) 


The  quantity  in  brackets  has  an  absolute  maximum  as  a function  of  P . 

Therefore,  above  a certain  value  of  j there  exists  no  steady-state  f or  T . 

sw  w 

At  this  critical  current  density,  the  temperature  and  electric  field  at  the 
cathode  continue  to  increase  until  arcing  occurs. 


From  the  form  of  Equation  (131)  it  is  clear  that  arcing  may  be  forestalled 
by  effective  cooling.  In  particular,  the  ratio  * /t  should  be  maximized  by 
using  thin,  highly  conductive  cathodes.  It  is  also  desirable  to  use  gases 
with  low  normal  current  densities  so  that  the  heat  loading  in  the  sheath  is 
minimized.  By  using  very  clean,  water-cooled,  metal  electrodes,  a normal 


glow  can  be  maintained  in  pure  hydrogen  at  pressures  as  high  as  13  atm 
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and  currents  as  high  as  14  amps. 

The  cathode  fall  potential  as  a function  of  reduced  current  density  in  nitrogen 
is  shown  in  Figure  25  for  a number  of  pressures.  The  cathode  is  1 cm  thick 
aluminum,  cooled  on  the  back  face  to  300°K.  The  curves  are  shifted  to  the 
left  at  high  pressures  because  of  the  reduced  gas  density  at  the  cathode  sur- 
face. At  the  critical  current  density  for  each  pressure,  the  slope  of  the 
curves  goes  to  infinity.  The  circles  indicate  the  points  where  the  wall  tem- 
perature reaches  the  melting  point  for  aluminum.  Obviously  the  emission 
process  will  become  thermionic  before  then  and  the  discharge  will  collapse 
into  an  arc.  At  pressures  above  one  atmosphere,  this  transition  occurs 
below  the  normal  current  density  and  the  discharge  always  arcs. 

8.  2 Field  Tailoring  and  Auxiliary  Ionization 

In  an  electric  discharge  with  transverse  gas  flow,  significant  heating  and 
chemical  reaction  may  occur  during  transit  through  the  device.  The  change 
in  gas  density  and  chemical  composition  will  alter  the  ratio  between  external 
and  volume  ionization  and  may  cause  instabilities  to  occur.  One  technique 
for  counteracting  this  effect  would  be  to  vary  the  electric  field  in  the  flow 
direction  by  segmenting  the  electrodes.  This  approach  would  also  allow 
individual  ballasting  or  current  control  of  each  segment  which  would  tend  to 
maintain  a uniform  current  distribution. 

The  changes  which  take  place  as  the  gas  flows  through  the  discharge  occur 
principally  in  the  positive  column.  Since  the  velocity  in  the  sheath  region 
is  small,  there  is  little  communication  between  upstream  and  downstream 
segments.  However,  the  gas  heating  which  occurs  in  the  sheath  causes  the 
thermal  boundary  layer  to  grow  at  a more  rapid  rate.  * The  expansion  of 
this  hot,  turbulent  gas  into  the  positive  column  together  with  the  normal 
bulk  heating  of  the  core  flow  leads  to  a severe  reduction  of  the  local  gas 
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density  and,  for  constant  field,  a corresponding  increase  in  E/N.  Most 
e -beam  sustained  discharges  are  operated  at  E/Ns  limited  by  some  form  of 
ionization  or  thermal  instability.  An  E/N  which  is  allowed  to  increase  in 
the  flow  direction,  then,  restricts  the  power  density  which  can  be  achieved 
in  these  devices.  If  the  electric  field  is  gradually  reduced  in  the  downstream 
direction,  then  theoretically  E/N  can  be  kept  constant  and  the  limiting  power 
density  increased. 

The  changes  which  occur  in  the  plasma  are  communicated  to  the  sheath  through 
the  local  current  density.  The  sheath  does  not  see  the  electric  field  in  the 
positive  column  or  the  gas  density  directly.  The  cathode  sheath  merely  has 
to  supply  the  current  required  by  the  discharge,  which  is  determined  by  the 
conductivity  in  the  positive  column.  If  the  local  conductivity  of  the  discharge 
drops  too  low,  then  the  current  density  required  from  the  cathode  may  exceed 
the  maximum  value  discussed  in  the  last  section  and  the  glow  will  collapse. 

On  the  other  hand,  there  are  instabilities  which  occur  in  the  positive  column 
which  lead  to  constriction  of  the  column  even  with  an  extended  cathode  glow'. 

The  presence  of  external  ionization  also  affects  the  cathode  sheath  only  in- 
directly. The  level  of  external  ionization  is  typically  small  compared  to 
volume  ionization  in  the  fall  region.  An  analysis  such  as  that  in  Section  7.2, 
including  an  external  source  term  in  the  electron  continuity  equation,  indicates 
that  the  electric  field  near  the  cathode  is  independent  of  the  source.  In  the 
transition  to  the  positive  column,  however,  the  electric  field  is  higher  than 
with  no  external  ionization.  This  high  field  in  a region  where  the  gas  density 
is  reduced  by  heating  does  strengthen  the  possibility  of  thermal  instabilities. 


CONCLUSIONS 


9.  0 

The  research  under  this  contract  has  answered  many  questions  relating  to 
basic  processes  in  the  plasma  sheath.  In  particular,  the  effects  of  nonequi- 
librium energy  distributions  and  gas  heating  on  sheath  characteristics  and 
stability  have  been  addressed.  However,  the  present  findings  do  not  close 
the  book  on  the  study  of  interactions  between  sheath  and  plasma  in  high  pres- 
sure discharges  and  ionized  gas  flows. 

Two  particular  aspects  of  the  sheath  problem  require  further  study.  One  is 
the  treatment  of  the  cathode  surface.  The  potential  drop  across  the  cathode 
sheath,  and  hence  the  heat  loading  in  this  region,  can  be  reduced  by  increas- 
ing the  effective  electron  yield.  This  can  be  accomplished  by  the  use  of 
oxides  as  dielectric  layers  or  the  use  of  helium  as  a gas  layer  above 
the  surface.  Furthermore,  the  tendency  of  the  positive  column  to  constrict 
and  eventually  collapse  at  high  pressures  may  be  averted  by  employing  re- 
strictive or  semiconducting  cathode  materials  to  provide  localized  ballasting. 
The  resistive  cathode  would  also  allow  tailoring  the  field  without  alternating 
metallic  and  dielectric  strips  which  result  in  strong  transverse  fields  at  the 
surface. 

The  other  area  needing  more  study  is  the  sheath -plasma  interface.  Under 
the  present  contract  it  was  not  possible  to  obtain  a smooth,  self-consistent 
field  transition  in  this  region  using  the  nonequilibrium  code  and  the  equilib- 
rium model  is  questionable  in  the  operating  regime  where  stability  becomes 
a factor.  The  problem  can  be  solved  by  modifying  the  iteration  scheme 
used  to  obtain  a self-consistent  field.  This  would  allow  the  calculations  to 
be  extended  to  the  positive  column  where  the  field  is  uniform  and  boundary 
conditions  can  be  established  without  the  presence  of  an  anode. 
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APPENDIX 


Asymptotic  Solutions 

The  collisional  Boltzmann  equation  for  electrons  suffering  elastic  and  in- 
elastic collisions  is  not  amenable  to  exact,  closed-form  solution  except 
in  the  case  of  constant  collision  frequency.  In  the  general  case,  numerical 
methods  are  employed  using  tabulated  experimental  cross  sections.  In  order 
to  make  the  solution  tractable,  it  has  usually  been  necessary  to  make  some 
approximation  in  the  angular  dependence  of  the  velocity  distribution  function. 
The  most  common  technique  is  to  expand  this  function  in  a series  of  Legendre 
polynomials,  which  is  then  truncated  after  two  terms.  This  approximation, 
which  has  been  labelled  PI,  is  very  accurate  at  moderate  electric  field 
strengths  in  cases  where  the  inelastic  cross  sections  are  not  too  large  and 
the  elastic  cross  section  is  isotropic.  In  the  present  work,  we  have  used 
an  approximation  where  electrons  moving  with  the  field  are  placed  in  one 
group  and  electrons  moving  against  the  field  in  another.  This  technique . which 
will  be  called  G2,  is  most  accurate  at  high  field  strengths  where  high  energy 
electrons  are  scattered  predominantly  in  the  forward  or  backward  directions. 


In  order  to  examine  quantitatively  the  accuracy  of  these  and  other  angular 

approximations,  we  now  consider  a situation  where  the  Boltzmann  equation 

has  a closed  form  solution  in  the  limit  of  high  electron  energy.  The  cross  sections 

are  assumed  to  be  constant  and  isotropic.  The  total  inelastic  cross  section 

is  denoted  by  Q and  the  total  elastic  cross  section  by  Q . The  Boltzmann 
h e 

equation  can  then  be  written, 


N ~ C°S  9 (“ll  + -f?)  = -Qf  + foQe 


e + e. 


f ((  + € ) Q (f  + £ ) ( ! ) 

o h h h 


where  Q - and  the  other  symbols  are  defined  as  in  Section  4.  1. 

Along  the  characteristic, where  + § is  equal  to  the  maximum  energy  of 
electrons  emitted  by  the  cathode,  there  are  no  in-scattered  electrons  and  the 
last  term  in  Eq.  ( 1 ) vanishes.  (This  term  is  also  negligible  in  the  homo- 
geneous case  when  E/N  is  small  or  Q is  large.)  Considering  f as  an  explicit 
function  of  5 alone,  we  can  write  Eq,  ( 1 ) as, 
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E df 

e — cos  0 — = 
N d? 


-Qf  + £ Q 
o e 


The  partial  derivative  in  (2)  can  be  expressed  in  the  coordinate  system 
(£  , x)  where  f - §/ cos  0 is  the  total  energy  and  x = cos  0,  i.  e. 


0 df  _ df  1 - x df 
cos  0 — - x — + 

a?  d€  2€  dx 


At  high  energies  the  second  term  on  the  right  is  negligible  compared  to  the 
first  except  in  a small  angular  range  around  x = 0.  Therefore,  we  can  write 
Eq.  (2)  in  the  asymptotic  limit  as, 


df 

x — = -f  + 
du 


2Q  J , 


where  u = NQc/eE. 


The  solution  of  this  equation  has  the  form 


f(u,  x)  = g(u)  h(x) 


where  the  functions  g and  h satisfy 


_L 


g 


2Q  J , 


hdx  = (1  - a x)h 


The  constant  a must  be  positive  to  meet  the  boundary  condition,  f -*  0, 
at  infinity. 


The  solution  of  Eq.  (5)  is 


g(u)  = e' 


h(x)  = 


Q / Q 

e 


1 - ax 


with  the  condition  that  — I hdx  = 1.  This  last  equality  fixes  a so  that 

2 *'-1 

there  is  a single,  well-defined  function  describing  the  electron  energy  dis- 


tribution, 


Q /Q 
e -au 

f(u,x)  = : — e 

1 - ax 


with  a given  by  a = —r  In  ( j + a 1 . Thus,  in  the  asymptotic  limit,  the 

20  ' 1 " a ' eE/N 

distribution  approaches  a Maxwellian  form  with  temperature  T = — — r — . 

6 '-*•  W 

The  angular  dependence  exhibits  a peak  in  the  forward  direction  which  becomes 
sharper  as  a -»  1.  The  dependence  of  a on  Q^/Q^  *s  given  in  Table  I. 

Table  I 


Q./Q 

h e 


f./f 
1 o 


Ut 

2  o 


Vf 

3  o 


iJi 

4  o 


IJt 

5 o 


0, 5030 


0,  542 


0.  195 


0.  060 


0,  3 

1.0 

10. 

0. 7521 

0.958 

0.999 

0.921 

1.  566 

2.  727 

0.  561 

1.  586 

4.  318 

0.  305 

1.427 

5.  834 

0.  156 

1.211 

7.295 

0.073 

0.987 

8.  717 

Another  method  of  solving  Eq.  (4)  would  be  to  multiply  by  successive  Legendre 
polynomials  and  integrate  over  x to  obtain  the  following  series  of  equations. 


t 


1 df1  Q 
3 du  o Q o 

? 9f,  Sf 

— + -2-  + f = 

5 du  du  1 


0 


(7a) 


(7b) 


3 

7 


= 0 . . . . 


(7c) 


where  the  functions  f , f , , f , . . . are  defined  by  the  expansion 

° 1 2 

f(u, x)  = £ f (u)  P (x).  If  the  series  is  approximated  by  the  first  two  terms, 
£ = 0 * 1 

i.  e. , fp  = 0 for  f 2 2.  Then  the  solution  of  ( 7 ) is  given  by, 


f 

o 


-0.  u 
1 


and  f 


1 


0.  e 

1 


-v 


where  a,  = y/ 3 Q / Q . 

1 h 

This  is  the  PI  approximation  commonly  used  to  solve  the  Boltzmann  equation 

in  the  general  case  of  nonconstant  cross  sections.  Once  again  the  solution  is 

c E / N 

Maxwellian  with  temperature  T^  = — ^ q the  value  of  CLj  approaches  the  exact 

result,  a,  at  small  values  of  Q^/Q^  the  limit  of  large  inelastic  cross 
section,  a goes  to  V~3  while  a goes  to  1.  Thus,  the  temperature  in  the 
tail  of  the  distribution  may  be  off  by  as  much  as  a factor  of  y/~T  . 
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If  three  terms  are  retained  in  the  angular  expansion  of  f,  the  solution  of  ( 7 ) 
is  given  by 


f = e 
o 


"a2U 


fl  = 


'a2U 


1 - 


4 2 

7?  °2 


a 1 


f2  = T 


-a2u 


1 ,-i-  a2 
15  2 


where  Ct^  - 


3 Q /Q 
n 


1 ♦ — Qh/Q 


This  is  the  P2  approximation. 


Note  that  a approaches  V 5/3  for  large  Q /Q  . 

2 he 


In  order  to  see  how  well  these  approximations  represent  the  exact  distribution, 

we  have  expanded  the  function  in  Eq.  (6)  in  a Legendre  series.  The  coefficients 

of  this  series  are  given  in  Table  I for  a number  of  values  of  Q /Q  . The 

h e 

coefficients  of  the  truncated  expansions  are  given  in  Table  II  along  with  the 

appropriate  cl  values.  As  the  inelastic  cross  section  increases  relative  to  Q^, 

more  terms  are  required  in  the  angular  expansion  to  adequately  represent  the 

distribution  function.  This  is  because  the  distribution  is  becoming  more  and 

more  peaked  in  the  forward  direction.  In  the  limit  as  Q -*  the  angular 

h 

distribution  approaches  a delta  function  which  has  the  Legendre  coefficients 
1,  3,  5 2C  + 1,  . . . 


Table  II 


Qh/Qe 

0.  1 

0.3 

1.0 

ai,2 

PI 

P2 

PI 

P2 

PI 

P2 

0.  522 

0.  516 

0.  832 

0.  764 

1.  224 

1.035 

Vfo 

0.  522 

0.  541 

0.  832 

0.905 

1.  224 

1.  449 

f2/fo 

0.  182 

0.  461 

1.  000 
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Before  leaving  the  angular  expansion  technique,  it  is  interesting  to  compare 
the  infinite  series  solution  with  the  exact  solution  obtained  by  separation  of 
variables.  By  a process  of  successive  substitution,  it  can  be  shown  that  the 
equations  ( 7 ) are  satisfied  by  functions  f^  which  have  a common  exponential 
factor  exp(-0ioou),  where  cx^.  is  given  by  the  continued  fraction, 


Q 
e 

Q 


a. 


1 - — 


2 

5 a«> 


2/3  a 


3/5  a. 


1 - — — a 

7 « 


4/7  a 


1 _ — aM  — 
9 00  1 


In  this  expression  the  two  series  are 


and 


n 

2n  + 1 


n 

2n  - 1 


The  value  of  found  in  this  way  is  equal  to  the  exact  value,  a , given  pre- 
viously. Therefore,  the  solution  using  the  Legendre  expansion  approaches  the 
exact  asymptotic  solution  as  more  and  more  terms  are  included. 

Another  very  simplified  approach  in  dealing  with  angular  velocity  distributions 
is  to  divide  the  electrons  into  two  groups.  One  moving  with  the  electric  field 
and  the  other  moving  against  it.  This  type  of  a distribution  is  most  nearly 
realized  in  the  cold  cathode  region  of  a glow  discharge  where  the  scattering 
of  high  energy  electrons  is  predominantly  forward  or  backward  and  the  inelastic 
cross  sections  are  comparable  to  the  elastic.  This  is  also  precisely  the  regime 
in  which  the  Legendre  expansion  has  the  most  trouble  converging.  In  this  case, 
the  Boltzmann  equation  is  written  as  follows, 


eE 

N 

+ . 
Sc 

-(q;  + <V£+  + Q;(- 

(8a) 

eE 

N 

Sf_ 

SC 

- + <q;  + Qh,f- 

(8b) 

where  Q“  is  the  cross  section  for  elastic  backs catte ring.  For  equal  forward 
and  backward  scattering,  Q~  = CWZ.  Notice  that  the  cross  section  for  for- 
ward scattering  does  not  appear  at  all  in  the  equations.  Only  the  backs  catte  ring 
and  inelastic  collisions  are  significant  in  shaping  the  distribution  function. 

The  solution  to  Eqs.  (8)  in  the  asymptotic  limit  is 


/ °h  -Bu 

f+  " f“  V 2Qe  + °h  6 

where  B = 70^20”  + Qh)  /Q  and  u = NQc/eE  as  before.  Once  again  the 
distribution  approaches  a Maxwellian  at  high  energy.  This  time  the  coef- 
ficient of  u in  the  exponent  goes  to  1 as  Q^/CT  increases  just  as  in  the  exact 
asymptotic  solution.  However,  for  small  values  of  Qh/Q~,  the  approximation 
is  not  as  good  as  PI.  This  is  because  of  the  conflicting  assumptions  of  isotropic 
scattering  in  the  one  case  and  forward-backward  scattering  in  the  other.  For 
real  gases  at  energies  above  10  eV,  the  cross  sections  are  more  nearly  des- 
cribed by  the  latter  assumption. 

It  would  be  useful  to  establish  a criterion  for  terminating  the  angular  expansion 
which  would  result  in  a good  approximation  to  the  distribution  function.  If  we 
require  that  the  term  which  is  dropped  in  the  series  of  Eqs.  (7)  be  less  than 
10%  of  those  remaining,  then  the  criterion  for  using  the  Pn  approximation  is 

fn  + 1 < o 1 n(Zn  + 3i 

f (n  + 1 )(2n  - 1) 

n - i 


JL 
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Using  Table  I,  we  see  that  PI  is  a good  approximation  for  Q^/Qe  5 0*  1 while 

at  least  seven  terms  must  be  retained  for  Q /Q  = 0.  3.  With  larger  inelastic 

h e 

cross  sections  the  criterion  is  never  satisfied. 


An  alternative  condition  is  to  require  that  f.  /f  in  the  truncated  expansion 

l o 

be  within  5%  of  the  exact  asymptotic  value.  This  criterion  is  satisfied  for 

PI  at  Q /Q  = 0.  1,  for  P2  at  0.  3,  and  for  P3  at  1.  The  error  in  a which 
he 

goes  as  one  over  the  electron  temperature,  is  about  equal  to  the  error  in 

f /f  , which  is  proportional  to  the  drift  velocity  of  the  high  energy  electrons. 

1 o 

Thus,  the  electron  transport  coefficients  are  much  more  accurately  predicted 
by  the  expansion  method  than  are  the  angular  distributions  themselves. 


GLOSSARY 


Subscript  denoting  attachment 


Attachment  Rate 


Electron  acceleration  vector 


Coefficients  for  Legendre  expansion 
Magnetic  field  vector 

Electron  velocity  vector  and  magnitude 


Heat  capacity  per  unit  mass 
y-component  of  electron  velocity 
Molecular  velocity  vector  and  magnitude 


Detachment  rate 


Anode -cathode  separation 


Cathode  fall  thickness 


Normal  cathode  fall  thickness 


Absolute  value  of  electronic  charge 


E,  E(y)  Electric  field  vector  and  magnitude 


Electric  field  at  cathode  surface 


f+'f- 


Electron  distribution  function 

Distribution  of  electrons  moving  with  and  against  the  electric  field 
Current  distribution  at  anode:  f = f(e,y  - d)  e 
Isotropic  distribution  function 

1/2 

Probability  distribution  at  anode:  F = f(e,y  - d)  e 


Magnitude  of  relative  velocity 
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Subscript  denoting  inelastic  collisions 


< 


i 


Je’J  + 


Jn 

k 

K 

K 

e 

m 


m+,  m_ 


M 

M 


Stagnation  enthalpy 
Subscript  denoting  ionization 
Radiation  intensity  at  cathode 
Total  current  density 
Electron  and  ion  current  density 
Normal  current  density 
Boltzmann's  constant 
In- scattering  operator 

In-scattering  operator  for  elastic  collisions 

Electron  mass 

Ion  masses 

Mach  number 

Molecular  mass 


n 


N 

N" 


P 

P„ 


q 

q. 


Electron  and  ion  number  densities 
Neutral  number  density 
Neutral  excited  state  density 
Pressure 
Reduced  pressure 
Discharge  power  density 
Legendre  polynomials 
Heat  flux 

Heat  flux  at  edge  of  cathode 

l 


1 1 


4 


Q Total  collision  cross  section 

Q Rate  of  heat  addition 

Q Total  attachment  cross  section 

a 

Q Total  elastic  cross  section 

e 

Total  inelastic  cross  section 

Q Total  ionization  cross  section 

i 

Q Total  recombination  cross  section 

r 

r Subscript  denoting  recombination 

r*  Position  vector 

r Electron-ion  recombination  coefficient 

e 

r Ion-ion  recombination  coefficient 

1 

s Subscript  denoting  superelastic  collisions 

S External  ionization  source  term 

t Time 

t Thickness  of  cathode 

T Gas  temperature 

T^  Adiabatic  wall  temperature 

T Temperature  of  back  face  of  cathode 

T Temperature  of  front  face  of  cathode 

w 

T Free  stream  temperature 

u,  v Flow  velocities  in  x,  y directions 

U Free  stream  velocity 
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Ion  drift  velocities 


vv- 


n 


W 


m 


Electron  drift  velocity 

Cathode  fall  potential 

Normal  cathode  fall  potential 

Electron  drift  velocity 

Position  coordinate  in  flow  direction 

Position  coordinate  parallel  to  electric  field 

Spherical  harmonics 

Ionization  rate 


i 

I 


a 


Y 


t 

l 

fi 

e” 


£ . 
1 

£ 

O 

3 


Townsend  ionization  coefficient 
Polar  scattering  angle 

Secondary  electron  emission  coefficient 

Photoelectric  efficiency 

Boundary  layer  thickness 

Boundary  layer  displacement  thickness 

Thickness  of  laminar  sublayer 

Total  kinetic  energy 

Electron  mean  energy 

Energy  lost  in  inelastic  collision 

Energy  lost  in  ionizing  collision 

Permittivity  of  free  space 

Angle  between  electron  velocity  and  -E 
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J 


i 


Thermal  conductivity  of  cathode 
Thermal  conductivity  of  gas 
Coefficient  of  viscosity 
Electron  mobility 
Kinematic  viscosity 
Ion  collision  frequency 
Kinetic  energy  in  field  direction 
Gas  Density 

Free  stream  gas  density 
Differential  elastic  cross  section 
Differential  inelastic  cross  section 
Recombination  cross  section 
Shear  stress 

Characteristic  times  for  gas  heating  in  the  cathode  sheath 

Electron  potential  energy  or  field  potential 

Work  function  of  cathode  material 

Azimuthal  angle  for  electron  velocity 

Viscous  dissipation  function 

Azimuthal  scattering  angle 

Solid  angle  for  scattering 
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